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ABSTRACT 


Interpolation  operations  have  wide  applications  in  signal  processing  systems. 

An  efficient  interpolation  technique  is  developed  and  demonstrated.  This  is  an  iterative 
technique  in  which  interpolation  output  at  any  stage  is  combined  with  the  input  to  that  stage 
so  that  the  interleaved  result  forms  the  input  to  the  succeeding  stage.  Sampling  rate 
increases  of  high  order  are  not  difficult  to  achieve.  Resampling  of  high-order  interpolator 
outputs  permits  precise  control  of  sampling  rate  conversion  ratios. 

From  this  basic  work  in  interpolation,  a digital  filter  synthesis  procedure  is  developed. 
The  resulting  class  of  finite  impulse  response  filters  is  competitive  in  speed  and  storage 
with  conventional  designs.  These  filters  have  the  property  that  the  passband  frequency 
function  of  a generic  design  can  be  scaled  by  any  desired  factor  in  frequency  while  main- 
taining stopband  suppression  to  specified  levels.  This  synthesis  procedure  uses  a 
"stretch-and-fill"  iteration  operation.  The  Impulse  response  into  a synthesis  stage  is 
stretched , and  the  vacated  spaces  are  then  filled  from  new  data  synthesized  at  that  stage. 

The  "stretch"  performs  the  frequency  scaling  function  while  the  "fill"  operation  eliminates 
image  responses  that  would  otherwise  appear.  Computation  of  the  coefficients  for  these 
filters,  of  any  order,  is  trivial. 

A FORTRAN  computer  program  which  uses  the  developed  algorithm  is  provided 
for  general  interpolation  service.  Auxiliary  use  of  this  program  for  obtaining  derived 
filter  coefficients  is  described,  and  examples  of  filter  synthesis  are  given.  Special  entry 
points  are  provided  for  data  format  mode  conversion:  complex-to-real  and  real-to-complex. 

Equivalent  procedures  for  sampling  rate  reduction  are  also  presented  and  discussed. 
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SECTION  I 


'NTRODUCTION 

The  basic  technique  of  Interpolation  Is  familiar  to  most  engineers  in  the  context  of 
"reading  between"  tabulated  points  in  tables  of  functions  or  other  data.  This  may  be  an 
unfortunate  introduction  to  a useful,  fascinating,  and  often  overlooked  technique  for  reasons 
that  go  beyond  the  immediate  mental  anguish  of  doing  the  proportions  correctly. 

Satisfactoly  linear  interpolation  requires  a fineness  in  the  spacing  of  the  abscissa 
values  or  "sampling  rate"  that  is  often  far  greater  than  that  required  to  contain  the  basic 
behavior  information  of  the  function  or  data.  This  is  especially  true  in  digital  signal 
processing  systems  where  band-limited  functions  having  accurately  known  spectra  are 
involved.  In  such  systems  economies  in  design  are  sometimes  overlooked  because  changes 
in  the  sampling  rate  within  the  system  would  be  required.  Such  changes  are  most  conven- 
iently  accomplished  by  the  resampling  of  interpolated  data.  Such  interpolation  must  be 
done  inexpensively  using  designs  which  allow  control  of  introduced  errors. 

In  some  cases  sampling  rates  that  are  adequate  to  permit  reconstruction  of  the  con- 
tinuous waveform  are  inadequate  for  certain  nonlinear  operations,  with  multiplication  being 
a common  example.  In  these  cases  interpolation  may  be  used  to  raise  the  sampling  rate 
just  prior  to  a nonlinear  operation  of  this  type. 

Interpolation  is  obviously  not  confined  to  the  time  domain.  In  the  spatial  domain  for 
example,  expensive  beamformer  processors  can  sometimes  be  saved  by  interpolation- 
synthesis  of  outputs  corresponding  to  intermediate  pointing  angles. 

Data  format  or  mode  conversion  between  "real"  data  and  "complex"  data  is  often 
required  in  signal  processing  work.  This  conversion  problem  can  be  shown  to  have  a simple 
interpolation  solution.  Media  conversions  often  involve  interpolation  because  of  differences 
in  appropriate  sampling  rates.  Samples  of  a waveform  taken  at  or  near  the  Nyquist  rate 
(two  times  the  highest  signal  frequency)  might  be  adequate  for  machine  processing.  However 
this  rate  would  be  totally  inadequate  for  display  use.  Try,  for  example,  to  mentally  re- 
construct a sine  wave  given,  say,  three-equispaced  randomly  located  sample  values  taken 
from  one  period.  Interpolation  prior  to  display  is  often  a necessity. 


1-1 


r i 

B > — 1 


In  this  paper  an  interpolation  technique  is  developed  which  is  b '■sic ally  quite  simple 
and  is  very  easy  to  use.  The  technique  appears  to  offer  computational  advantages  over  more 
conventional  interpolation  methods.  A modification  of  the  interpolation  process  leads  to  a 
filter  synthesis  procedure  which  is  described.  Applications  such  as  mode  transformations 
are  treated  in  detail.  A computer  program  for  interpolation  and  mode  conversion  is 
presented  along  with  instructions  for  use. 

The  basic  techniques  developed  for  interpolation  may  also  be  used  for  sampling-rate 
reduction.  Computer  programs  and  design  procedures  are  presented  which  permit  con- 
venient sampling  rate  conversions  either  up  or  down. 
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SECTION  n 

BASIC  SINGLE-LEVEL  INTERPOLATION 

2.1  FREQUENCY  DOMAIN  APPROACH  TO  INTERPOLATION 

For  present  purposes  it  is  convenient  to  formulate  the  interpolation  problem  in  the 
context  of  the  tapped  delay  line  representation  of  Figure  2-1.  A sinusoidal  signal 

s (t)  = cos  27r  ft  (2-1) 

is  applied  to  a delay  line  having  2N  taps  spaced  T-seconds  apart.  In  later  development 
T will  represent  a sampling  period  for  a sampling  rate  R where 


It  is  desired  that  the  output  at  the  center  of  the  line  be  obtained  where 

c(t)  = cos  2ir  f [ t - (2N  - 1)  -y]  (2-3) 

An  estimated  center  output  c (t)  is  to  be  delivered  from  the  weighted  sum  of  the  tap  voltages 
so  that 


N 

= S I aR  cos  2 v { [t  + (2k  - 1)  y-  J 

+ a_k  cos  2ir  f |t  - (2k  - 1)  -y  J j 


using  the  signal  at  the  center  tap  as  reference.  Expansion  and  rearrangement  of  the 
terms  in  Equation  (2-4)  yields 


■[& 

-ft 


(ak  + a_k)  cos  2tt  f (2k 


cos  2n  ft 


(afc  - a_k)  sin  2n { (2k  - 1)  j sin  2n  ft 
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a.  Tapped  Delay  Line  Representation 


-3 


b.  Phasor  Diagram 


Figure  2-1.  Interpolation  Processor,  Continuous  Time  Version 


Since  the  center  tap  signal  was  used  as  reference  for  Equations  (2-4)  and  (2-5),  we  may 
ignore  the  (2N- 1)  T/2  delay  term  of  Equation  (2-3)  and  represent  cos  2 r ft  as  the  desired 
output  signal.  It  then  becomes  apparent  that  the  first  term  of  Equation  (2-5)  represents 
an  in-phase  component  while  the  second  term  represents  a quadrature  component  relative 
to  the  desired  cos  2 ir  ft  signal.  In  order  to  eliminate  the  undesired  quadrature  term  choose: 


and  define 

\ ' 2ak 


So  that  Equation  (2-5)  becomes  after  Equations  (2-2),  (2-6)  and  (2-7)  are  used 


V°<M2k-l>  i- 


COS  2 IT  ft 


The  estimate  of  Equation  (2-8)  has  no  phase  error  so  that  the  bracketed  factor 
may  be  interpreted  as  the  real  filter  transfer  function  H (f)  of  the  estimator  of  Figure 
2-la  (time  delay  of  (2  N - 1)  T/2  understood)  where 


H (f)  = Lj  bk  cos  (2k  - 1)  y ,rT  (2_9) 

k = 1 [2  ) 

If  H (f)  were  equal  to  unity,  the  estimation  would  be  perfect.  As  shown  in  Figure  2- lb  the 
coefficient  choice  of  Equation  (2-6)  Insures  that  the  resultant  phasor  c lies  along  the  true 
center  phasor  line  c.  However,  the  amplitude  of  c may  differ  from  unity  and  hence  intro- 
duce estimation  errors.  The  problem  then  becomes  one  of  choosing  the  coefficients 
bk  so  that  H (f)  of  Equation  (2-9)  remains  as  close  as  possible  to  unity  over  a frequency 


range  which  will  be  confined  to 

0 ^ f < 


(2-10) 
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Both  the  phasor  diagram  of  Figure  2-lb  and  Equation  (2-9)  show  the  difficulty 
involved,  namely  that  the  component  angular  relationships  are  functions  of  frequency. 

At  zero  frequency  the  phasors  are  all  collinear,  and  the  sum  of  the  coefficients  should 
be  set  equal  to  unity.  As  the  frequency  is  increased  the  phase  relationships  change  and  a 
different  coefficient  set  is  required  for  unity  gain.  Since  the  single  frequency  analysis 
is  meant  to  represent  performance  at  only  one  of  a band  of  signal  frequencies,  the  choice 
cannot  be  changed  aa  a function  of  frequency.  An  optimum  bfc  set  (in  some  sense)  must  be 
found  for  the  signal  characteristics  of  the  application. 

2.2  PROCESSOR  COEFFICIENT  DETERMINATION 

Examination  of  Equation  (2-9)  reveals  that  an  odd  harmonic  series  of  cosine  terms 
is  Involved  with  each  cosine  term  having  unit  value  at  f = 0 and  zero  value  at  f = (R/2). 

Thus  H (f  = R/2)  = 0 regardless  of  the  bfe  set  chosen.  Define  now  a maximum  value  of 
f/(R/2)  and  consider  design  procedures  confined  to  the  range 


0 s 


s y 


(2-11) 


It  is  desired  to  minimize  the  peak  error  (deviation  from  unity)  of  H(f).  If  the  error  is 

defined  as  6 then  the  quantity  to  be  minimized  is  6 where 

max 

6 = |l  - H (f)|  <2“12> 

That  is,  we  seek  that  set  |b^  which  minimizes  the  maximum  value  of  6 of  Equation  (2-12) 
for  a given  y , N.  As  expected,  the  computational  load  (required  N value)  is  increased 
when  <5  max  is  made  smaller  and  when  y is  made  larger  (closer  to  unity). 

Consider  now  the  trivial  but  illustrative  case  in  which  N = 1.  H(f)  now  becomes 

H(f)  = bj  cos  j-  (2-13) 

T 

When  f = 0,  H(f)  is  equal  to  bj  and  when  f/(R/2)  = y , H(f)  » b^  cos  w y/2.  Since  the  H(f) 
behavior  is  monotonic,  b^  must  be  chosen  to  equalize  the  errors  at  the  frequency  extremes. 
The  choice  is  clearly 


bl  " 


1 + cos 


(2-14) 
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and 


_ 1 - cob  ffy/2 
0 max  1 + cos  * y/2 

which  may  be  rewritten  as 

,«  1-6  max 

008  ff?/2  = 1+6  max 


(2-15) 


(2-16) 


The  interaction  of  performance  and  bandwidth  for  this  simple  case  is  easily  calculated. 

For  6 levels  of  -20,  -40,  and  -60  dB,  the  respective  y values  are  0. 39,  0. 13,  0. 04. 
m&x 

Thus  for  a fixed  computation  load  (N  = 1 in  this  case)  one  can  only  trade  performance  for 
bandwidth. 

The  above  procedure  represents  a modified  form  of  linear  interpolation  using  data 
from  the  delay  line  tap  pair  adjacent  to  the  center  point  of  Figure  2-la.  Ordinary  linear 
interpolation  would  involve  (a^  a ) values  of  1/2  and  a unit  value  of  b^  Note  that  bj  from 
Equation  (2-14)  is  always  greater  than  one.  This  departure  from  ordinary  linear  interpolation 
is  necessitated  by  the  desire  to  minimize  the  peak  error  in  H (f)  over  the  operating  bandwidth. 
The  result  in  this  case  is  a 6 dB  reduction  in  peak  error  level  as  compared  to  the  normal 
linear  interpolation  usage  of  the  tap  signals. 

The  computational  approach  used  for  general  {b^}  determination  will  now  be  described. 
First  rewrite  Equation  (2-9)  letting 

x - (2-17) 


so  that 


= b.  cos  (2  k - 1)  x 


(2-18) 


In  the  x interval  from  0 to  y select  a set  of  equally  spaced  points  { x^}  and  consider  the 
weighted  squared  error  sum 


wi  [' -«<*,>] 


(2-19) 


The  process  starts  with  unit  weights  and  a {b^}  set  is  determined  by  the  usual  (Ref.  1) 
least  squares  procedures  In  which 

cos  (2  k - 1)  -y-  x 

la  used  to  define  the  coordinate  functions  of  the  approximation.  After  the  least-squares 
operation  is  complete,  the  weights  Wj  are  modified  In  proportion  to 

|l-H(x1)| 

so  that  the  regions  of  high  error  are  emphasized.  Then  a new  weighted  least-squares  fit 
procedure  is  Invoked  to  derive  a modified  coefficient  set  { b ^ } . After  a few  such  Itera- 
tions It  will  be  found  that  the  peak  errors  of  H (x)  over  the  approximation  region  are 
equalized  and,  coincidentally,  minimized.  (Several  years  ago  the  author  wrote  an  inter- 
active time-share  program  for  this  and  a variety  of  other  curve-fitting  problems  which 
leans  heavily  on  a Honeywell  H-635  auxiliary  library  routine  named  LSQMM.  The  names 
M.  A.  Martin  and  F.  E.  Liller  are  associated  with  this  program. ) This  procedure  appears 
to  have  more  general  applicability  and  may  be  found  more  convenient  to  use  than  the 
REMEZ  exchange  algorithm  employed  for  similar  purposes  (Ref.  2) 

These  numerical  procedures  were  used  to  obtain  the  performance  data  shown  in 
Figure  2-2.  The  particular  cases  displayed  here  were  run  for  specific  purposes  to  be 
described  in  later  sections.  For  the  present,  the  results  of  Figure  2-2  may  be  discussed 
in  terms  of  the  general  interpolator  design  problem. 

The  time-share  program  previously  mentioned  was  run  for  a variety  of  y values 

over  an  appropriate  set  of  N values.  For  each  N value  the  peak  error  In  H (x)  was  found 

and  expressed  in  decibels.  Plots  of  6 in  dB  as  a function  of  N were  then  made  using 

max 

6 as  a parameter.  For  a specific  y , say  80%,  the  peak  error  decreases  monotonically 
with  N as  expected.  For  this  particular  y only,  the  error  performance  for  least-squares 
determination  of  the  coefficient  set  {b^}  is  shown  circled  above  and  to  the  right  of  the 
normal  minmax  performance  curve  marked  "y  = 80% ".  It  appears  that  the  minmax 
iteration  is  worth  the  order  of  6 to  8 dB  in  additional  error  suppression. 

The  rapidly  escalating  processor  costs  incurred  as  y approaches  one  may  be  seen  by 
looking  along  the  -35  dB  error  line,  for  example.  At  y = 80%  an  N of  6 will  suffice, 
at  y = 90%  an  N of  11  is  indicated,  while  at  y = 95%  the  computation  price  has  risen  to 
N = 23.  For  the  smaller  y values  N goes  down  rapidly,  with  simple  two-point  interpolation 
(N  = 1)  providing  error  performance  in  excess  of  -60  dB  for  y of  values  of  about  15%  and 
less. 


In  the  discussion  up  to  this  point  the  Interpolation  problem  has  been  presented  in 
the  context  of  continuous  time  functions  made  available  only  at  equispaced  delay  intervals 
T.  Simple  modifications  of  the  structure  of  Figure  2-la  provide  a convenient  conversion 
to  the  discrete  time  case  of  interest  in  digital  signal  processing.  As  shown  in  Figure  2-3 
the  delay  line  input  signal  is  sampled  at  rate  R = 1/T  by  the  Impulse  train  p (t).  Additional 
changes  from  Figure  2-la  involve  summation  of  the  c (t)  output  with  the  output  of  the  a_^ 
tap  after  this  tap  output  pulse  is  delayed  by  T/2  seconds.  At  any  given  input  pulse  time, 
tQ,  the  output  of  tap  a_1  represents  input  signal  as  it  existed  (N  - 1)  T seconds  prior  to  tQ. 
The  output  c approximates  input  signal  as  it  existed  (N  - 1)  T + T/2  seconds  prior  to  t . 

A ^ 

Since  the  c sample  represents  an  earlier  epoch  of  s (t)  than  the  a_^  sample,  the  latter  must 
be  delayed  by  T/2  to  maintain  proper  order.  Thus  the  structure  of  Figure  2-3  may  be 
considered  to  operate  on  input  data  sampled  at  rate  R and  produce  at  its  output  sampled 
data  at  rate  2R.  Note  also  that  the  delay  through  the  processor  of  Figure  2-3  is  (N  - 1) 

T + T/2  seconds  or  2N  - 1 output  rate  sampling  periods. 


In  the  analysis  of  Figure  2-3  which  follows  it  is  assumed  that  s (t)  has  a power 
density  spectrum  S (f)  where 

S (f)  = 0 for  | f | 2R/2  (2-20) 

Since  the  sampling  period  is  precisely  equal  to  the  delay  line  tap  spacing,  all  line  output 
voltages  will  be  zero  except  for  times  which  are  multiples  of  the  sampling  period  T.  The 
tap  voltage  (t)  will  be 

dx  (t)  = c(t)  p (t)  (2-21) 

The  output  on  the  a . tap  represents  signal  data  that  is  from  a later  time  than  the  data  on 

^ ™ A 

the  c(t)  bus  by  T/2  seconds.  For  this  reason  the  a_1  impulse  samples  of  Figure  2-3  must 
be  delayed  as  shown  in  order  to  maintain  proper  output  sequence.  This  T/2  delay  of  the 
a_^  sequence  is  exactly  equivalent  to 


d2  (t)  = s (t  - (2N  - 1)  T/2)  p (t  - T/2) 

The  impulse  train  p (t)  may  be  represented  by  the  Fourier  series 

oO 

p(t)  = ^ cos  2ir  nR  t 

n = 0 

so  that 

00 

P(t  - T/2)  = £ (-1)“  cos  2ir  nRt. 

n = 0 

Equation  (2-22)  may  be  written  as 

00 

d„(t)  = s (t  - (2N  - 1)  T/2)  £ (-1)11  cos  2x 


(2-22) 


(2-23) 


(2-24) 


nRt 


(2-25) 


n = 0 


and  the  other  adder  input  will  be 


dj(t)  = c(t)  p(t)  = c(t)  cos  2ir  nRt. 


2-26) 
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A sketch  of  Equation  (2-30)  for  an  S (f)  of  unity  in  the  range  | f ] - R/2  is  shown  in 
Figure  2-4a.  The  sampled  data  spectrum  at  the  input  to  the  interpolator  has  a periodicity 
of  frequency  span  R,  the  input  sampling  rate.  The  spectral  periodicity  at  the  interpolator 
output  has  a span  2R  consistent  with  the  doubling  of  the  sampling  rate  as  explained  relative 
to  Figure  2-3.  Thus  the  new  foldover  frequency  (one-half  the  sampling  rate ) is  R. 


Relative  to  this  new  foldover  frequency  we  may  define  an  effective  frequency  function 
G (f)  for  the  interpolation  processor  of  Figure  2-3 


G (f)  = |l  + H (f)|  , 0 <f<R/2 

G (f)  = |l  - H (R  -f)|,  R/2  <f  SR 
remembering  that 


(2-31) 


H (f  = R/2)  = 0 (2-32) 

The  complication  introduced  by  the  doubling  of  the  sampling  rate  from  the  sampled 
data  input  of  Figure  2-3  to  the  output  at  d (t)  may  be  resolved  by  assuming  an  input  data 
rate  of  2R  with  alternate  sample  values  all  equal  to  zero  (Ref.  3).  This  construction  leaves 
the  physical  input  signal  to  the  line  of  Figure  2-3  unchanged,  yet  allows  input  and  output 
power  density  spectra  to  be  directly  compared  in  the  context  of  an  equivalent  filter  G (f) 
of  Equation  (2-31). 


As  discussed  earlier,  an  ideal  H (f)  function  would  have  unit  value  out  to  the  original 
foldover  frequency  R/2.  The  resulting  G (f)  would  have  a gain  of  2 from  0 to  R/2  and  be 
zero  elsewhere.  This  would  produce  block  spectra  in  Figure  2-4  of  width  R centered 
at  even  multiples  of  R.  This,  of  course,  is  the  spectrum  which  would  result  if  the  original 
s (t)  were  directly  sampled  at  the  rate  2R. 

The  filter  function  G (f)  of  Figure  2-4a  and  Equation  (2-31)  may  be  discussed  in 
terms  of  a passband  from  0 to  y(R/2),  a transition  region  from  y (R/2)  to  (2-y)(R/2),  and 
a stopband  from  (2  - y)  (R/2)  to  R.  The  passband  and  stopband  are  directly  related  since 
Equation  (2-31)  shows  that  a deviation  6 of  H (f)  from  unity  in  the  passband  at  frequency 
f results  in  a response  of  magnitude  1 6 1 in  the  stopband  at  frequency  R - f.  Interpolation 
errors  result  in  the  desired  signal  being  at  an  Incorrect  level  at  the  output  (passband) 
accompanied  by  an  image  signal  created  within  the  processor  (stopband). 
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The  gain  of  2 in  the  passband  indicates  a reinforcement  of  desired  output  signal  from 
the  interpolated  series  by  the  signal  from  the  "direct-through"  series  from  tap  a ^ of 
Figure  2-3.  Remember  that  the  error  performance  calculations  shown  in  Figure  2-2  were 
made  relative  to  the  interpolated  series  only.  Thus  the  relative  error  values  of  the 
sampled-data  version  of  Figure  2-3  are  reduced  by  6 dB  from  the  values  given  in  Figure  2-2. 

The  transition  region  of  Figure  2-4  is  obviously  a region  of  poor  performance.  The 
desired  signal  level  errors  can  be  large,  resulting  in  high  levels  of  image-signal  generation. 
Error  effects  resulting  from  this  region  can  be  reduced  either  by  making  y large  (which 
can  be  computationally  very  expensive)  or  by  keeping  the  signal  spectrum  in  the  region 
y R/2  to  R/2  low.  The  latter  condition  results  when  the  sampling  rate  is  increased 
beyond  twice  the  low-pass  filter  cutoff  frequency  when  the  Initial  analog-to-digital  conversion 
is  made,  for  example. 

Filters  of  this  type  are  also  discussed  in  reference  5 and  are  described  as  "half-band 
non-recursive"  filters. 
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SECTION  III 
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I 


MULTILEVEL  INTERPOLATION 

3.1  EXTENSION  OF  TECHNIQUE  TO  GENERAL  L-1.EVEL  CASE 

If  the  sampling  rate  is  to  be  increased  by  a factor  of  2^,  an  iteration  of  the  technique 
of  Figure  2-3  may  be  used  to  good  advantage  as  shown  in  Figure  3-1.  The  general  L-level 
procedure  will  be  outlined  first,  with  detailed  design  considerations  to  follow.  (Different 
design  approaches  to  multistage  interpolation  may  be  found  in  references  7,  8 and  10. ) 

Starting  at  the  top  of  Figure  3-1,  we  reason  that  the  interpolated  data  at  the  output  of 
level  1 represents  the  best  estimate  of  intermediate  data  values  obtainable  from  the  alloca- 
ted processing  investment  made  at  that  level.  Thus  it  seems  eminently  reasonable  to  use 
these  intermediate  results  as  input  for  the  next  processing  stage. 

The  design  of  the  next  stage  (level  2 in  this  case)  may  be  relaxed  as  compared  to  that 
of  level  1 because  the  relative  signal  bandwidth  has  been  cut  in  half  by  the  1:2  increase  in 
sampling  rate  at  the  previous  level. 

For  example,  consider  an  error  level  line  drawn  horizontally  in  Figure  2-2  through 

the  y = 60,  30,  15, % curves.  If  the  signal  bandwidth  requirement  calls  for  a y = 60% 

first  level,  a proper  lst-level  N value  may  be  chosen  by  noting  the  intersection  of  the  error 
line  with  the  curve  y = 60%. 

Now  at  the  second  level  the  effective  y is  reduced  to  30%  and  a lesser  N will  be  satis- 
factory which  means  a lighter  processing  load  as  compared  to  level  1.  The  next  step 
results  in  a level  3 y of  15%,  and  so  on.  Eventually  the  working-level  y is  reduced  to  the 
point  where  N = 1 which,  in  effect,  means  linear  interpolation.  At  this  point  (say  at  the 
output  of  level  K)  the  procedure  is  changed  and  the  final  processor  simply  does  linear 
interpolation  of  2 -1  points  between  each  input  pair.  This  procedure  is  efficient  in  that 

only  that  amount  of  processing  needed  for  error  control  is  used  at  each  level. 

The  procedure  may  obviously  be  extended  to  provide  arbitrary  sampling  rate  multi- 
plication at  the  linear  interpolation  output  level. 

The  problem  specifications  determine  how  many  Figure  2-3  type  processors  will  be 

L-K 

needed  before  reaching  the  2 linear  interpolation  processor.  Once  the  linear  inter- 
polation level  is  reached  sampling  rate  multiples  become  inexpensive. 

The  above  discussion  is  conceptually  correct,  but  important  design  details  which 
were  passed  over  will  be  considered  in  the  next  two  sections. 
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Figure  3-2.  System  Filter  Functions  and  Spectra,  Critically-Sampled  Case 


3.2  DESIGN  CONSIDERATIONS.  CRITICALLY  SAMPLED  CASE 

The  output  spectrum  of  Figure  2-4a  was  based  on  a uniform  Input  signal  spectrum 

extending  to  R/2.  Deviations  of  this  curve  from  a repeated  block  spectrum  on  2R  Hz 

centers  represents  error  power.  A normalized  error  spectrum  Is  sketched  In  Figure  2-4b. 
o 

The  peak  level  6 Is  reached  In  the  pass  and  stopbands.  In  the  transition  region  a high 
density  level  is  always  reached.  Error  contributions  in  this  region  can  only  be  reduced 
by  narrowing  the  transition  region  (Increasing  y ) or  by  lowering  the  spectral  level  of 
the  signal  in  the  foldover-frequency  region  (oversampling). 


Error  control  In  the  multilevel  interpolator  Is  strongly  influenced  by  the  specific 
application.  In  much  of  the  previous  literature  (Ref.  3)  Interpolation  is  done  by  passing  a 
digital  signal  through  a digital  low-pass  filter  which  Is  characterized  by  a passband,  a 
transition  region  and  a stopband.  In  order  to  compare  Interpolation  methods,  the  design 
procedures  in  this  section  will  also  be  based  on  an  extended  stopband  in  which  uniformly 
low  error  level  is  maintained.  In  the  multilevel  interpolator,  error  spectrum  distribution 
can  take  many  forms.  One  form  arises  from  the  use  of  a low-pass  filter  in  the  operation. 

It  is  not  obvious  that  this  is  necessarily  the  best  error  distribution  for  all  applications. 

Figure  3-2(a)  shows  a uniform  input  signal  spectrum  extending  to  R/2.  Figure  3-2(b) 
shows  the  first-level  filter  function  and  output  spectrum  for  the  critically  sampled  case.  A 
y ^ is  chosen  at  the  first  level  which  yields  a transition  region  from  y ^ (R/2)  to  (2  - y1)R/2 
at  first-level  output.  The  error  spectrum  is  not  shown.  However,  it  is  known  to  peak  at 
the  center  of  each  transition  region.  This  first-level  design  will  produce  stopbands  in  the 
spectrum  as  shown  by  the  circled  number  ones,  indicating  stop  regions  contributed  by 
level  1. 
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The  level  2 design  uses  an  "extended"  y2  of  (2  - y^/2  in  order  that  the  created 
stopbands  shown  as  circled  two's  encompass  the  transition  region  and  error-spectrum 
peaks  associated  with  the  spectra  centered  at  multiples  of  2R.  Figure  3-2(c)  shows  the  stop- 
bands  produced  by  the  level  2 processor  while  Figure  3-2(d)  shows  the  level-two  output  spec- 
trum with  the  stop  regions  identified  as  to  the  processor  level  responsible. 


Starting  with  the  third  level,  the  y values  can  be  simply  one-half  the  previous  level 
value  as  shown  in  Figure  3-2(e)  (since  the  transition  regions  of  these  processors  will  corres- 
pond to  previously  produced  null  regions).  The  spectrum  and  null  region  identification  after 
three  levels  is  shown  in  Figure  3- 2(f).  Using  this  design  approach,  the  error  in  the  final  out- 
put is  essentially  the  error  contributed  by  the  first  level  processor.  The  first-level  error 
spectrum  peak  Is  carefully  maintained  (unfortunately)  in  order  that  all  other  error  peaks 
in  the  null  region  will  be  suppressed. 
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3.3  DESIGN  CONSIDERATIONS.  OVERSAMPLED  CASE 


A far  more  pleasant  design  task  is  presented  when  oversampling  is  done  so  that  a 
null  region  in  the  input  spectrum  exists  as  shown  in  Figure  3-3(a).  Now  the  value  can 
be  set  near  the  signal  cutoff  frequency  and  the  y value  can  be  halved  between  all  levels. 

The  transition  region  of  the  first-level  processor  falls  in  a null  region  created  by  the  signal 
which  is  identified  with  a circled  S.  The  remaining  plots  in  Figure  3-3  are  self-explanatory 
based  on  the  Figure  3-2  discussion.  There  are  no  error  spectrum  peaks  which  need  special 
attention  in  the  oversampled  case. 

If  the  signal  spectrum  does  not  extend  to  zero  frequency  further  design  economies 
are  possible.  In  the  coefficient  determination  process  described  in  Section  n,  par.  2.  2, 
me  region  over  which  H (f)  is  to  be  held  to  unity  is  now  confined  with  "don't  care"  regions 
both  above  and  below  the  passband.  This  results  in  better  passband  performance  for  a given 
N value,  and  corresponds  to  the  "stopband"  filtering  procedures  discussed  by  Schafer  and 
Rablner  in  Reference  3. 

In  the  discussion  relative  to  Figures  3-2  and  3-3  emphasis  was  placed  on  the  stop- 
bands  contributed  at  each  level.  The  transition  region  responses  have  cumulative  filtering 
effects  which  are  significant  as  will  be  noted  later  in  the  discussion  of  filter  synthesis. 


SECTION  IV 

FILTER  SYNTHESIS  APPLICATIONS 

4.1  PRELIMINARY  CONSIDERATIONS 

When  the  interpolation  process  discussed  in  Sections  n and  III  is  examined  as  a digital 
processor  and  the  equivalent  frequency  function  is  studied,  some  interesting  results  occur. 

The  work  in  this  section  leads  to  a class  of  digital  filters  having  several  interesting  properties: 

• The  impulse  response  has  a central  term  of  unity  and  has  a symmetric 
distribution  about  this  central  value.  The  resulting  finite  impulse  response 
(FIR)  filter  provides  a linear  phase  delay. 

• The  resulting  passband  region  frequency  function  can  be  scaled  arbitrarily 
by  stretching  the  given  impulse  response  and  filling  in  the  vacated  terms 
with  new  values. 

• This  scaling  operation  can  be  extended  without  limit  in  the  synthesis  procedure. 

The  numerical  computation  involved  in  coefficient  determination  is  trivial. 

• A large  variety  of  filter  designs  may  be  obtained  from  relatively  few  stored 
"seed"  parameters. 

• The  designer,  in  addition  to  scaling,  has  control  of  passband  width,  transition 
region  width  and  stopband  peak  levels.  Passband  ripple  is  not  directly  controlla- 
ble but  is  related  in  a complex  way  to  stopband  attenuation.  For  the  usual 
stopband  levels,  passband  ripple  is  not  a problem. 

• Simple  procedures  exist  for  controlling  passband  ripple  in  those  special 
cases  in  which  passband  gain  must  be  precisely  controlled. 

• Speed  and  storage  needs  of  these  filters  appear  to  be  competitive  with  the 
requirements  of  more  conventional  filter  designs. 

A brief  review  of  the  details  of  the  interpolation  operation  will  serve  to  Introduce 
the  filter  synthesis  work.  The  single  level  Interpolator  of  Figure  2-3  creates  two  inter- 
leaved outputs.  The  original  input  signal  (delayed)  from  the  a_^  tap  and  the  interpolation 
signal  on  the  d1  bus.  If  this  processor  is  cascaded  with  a similar  one  as  per  Figure  3-1, 
the  output  of  the  second  processor  will  contain,  intact,  the  output  of  the  first  processor 
interleaved  with  new  data  generated  by  the  second  processor.  The  output  of  the  second 
processor  will,  of  course,  also  contain,  intact,  the  input  to  the  first  processor  as  is 
obvious  by  induction.  The  data-rate  doubling  at  each  stage  makes  room  for  the  old  (input) 
to  be  interleaved  with  the  new  data  generated  at  that  stage. 
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In  the  general  case  of  an  L-level  processor  as  per  Figure  3-1  one  can  easily  show 
that  the  delay  through  the  processor  in  terms  of  output  pulses  is  simply 

L 

Dt  = £ 2 (L  + 1 ‘ k)Nk  - <2L  - 1)  = Nf  (4-1) 

k = 1 

where  Nk  is  the  number  of  coefficients  at  level  k (one-half  the  processor  length).  For 
L > 5,  the  linear  interpolation  operation  of  program  C432  always  inserts  a leading  zero  into 
the  impulse  response  of  the  equivalent  interpolation  filter.  This  produces  a delay  which  is 
one  greater  than  that  given  by  Equation  (4-1).  This  result  is  predicated  on  the  use  of 
= 1 for  k > 5 when  the  calculation  is  performed.  The  number  of  coefficients,  Nf,  on 
each  side  of  the  central  (unity)  term  of  the  impulse  response  of  the  interpolation  cascade  is 
also  given  by  the  Equation  (4-1)  calculation.  This  will  be  discussed  in  more  detail  at  a 
later  point.  The  unit  correction  for  for  L > 5 does  not  apply  to  the  Nj  calculation. 

In  terms  of  processor  load  it  is  convenient  to  define  an  "operation"  as  two  additions 
plus  one  multiplication.  Since  am  = a_m  in  these  designs,  one  sums  pairs  of  tap  voltages, 
multiplies  this  sum  by  the  coefficient,  then  sums  to  the  bus.  It  is  easy  to  show  that  the 
total  number  of  interpolator  operations  per  input  sample  is 

L 

O = £ 2(k_1)N.  (4-2) 

k = 1 

Consider  now  a level  one  processor  of  the  Figure  2-3  type  operating  in  the  context 
of  Figure  3-1  with  a ylt  of,  say,  80%.  The  spectrum  of  the  level  one  output  would  appear  as 
shown  in  Figure  3-2(b)  where  the  sampling  rate  is  2R  as  compared  to  R at  level  one  input. 
The  spectrum  of  Figure  3-2(b)  also  represents  the  frequency  function  (Figure  2-4a)  of  the 
level  one  processor  in  the  usual  "Fourier  transform  of  the  impulse  response"  context, 
based  on  the  output  sampling  rate  2R.  Note  that  as  a digital  filter,  the  level  one  processor 
has  a passband  edge  which  is  40%  of  foldover  or  one-half  the  input  design  value. 

The  input  sampling  rate  at  each  level  becomes  the  foldover  frequency  at  the  output 
of  that  level.  Thus,  the  second  level  filter  sees  an  entire  frequency-pattern  period  at  its 
input  transformed  to  fall  within  its  foldover  interval  as  defined  at  the  output  sampling  rate. 

In  Figure  3- 2(c)  an  extended  y is  used  to  place  the  entire  spectral  group  centered  at  2R  in 

i — — 
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In  Figure  3-2  the  frequencies  shown  are  absolute,  but  the  y relationships  to  foldover 
are  referenced  to  input  sampling  rate  values  as  was  appropriate  for  the  original  interpolation 
discussions.  For  present  purposes  reference  to  output  sampling  rates  is  more  appropriate, 
which  involves  simply  a scaling  factor  of  one-half.  At  the  output  of  the  third  level,  for 
example,  the  absolute  passband  edge  frequency  is  still  y ^ (R/2)  but  the  sampling  rate  here 
is  8R  which  means  that  the  effective  y has  been  cut  by  a factor  of  eight.  The  passband  edge 
at  the  output  of  level  three  is  only  10%  of  foldover  relative  to  the  level  three  output  sampling 
rate. 

Comparing  the  original  input  signal  spectrum  of  Figure  3-2(a)  with  the  level  three  out- 
put spectrum  of  Figure  3-2(f)  suggests  a digital  filter  with  a frequency  function  in  the 
Figure  2-4a  context  with  a y 1/8  passband  relative  to  foldover.  The  suggested  low-pass 
filter  function  can  be  rationalized  relative  to  the  sampling  rate  change  by  assuming  an  8R 
input  sampling  rate  with  seven  zero  samples  following  every  actual  sample.  The  extended 
y choice  at  level  two  of  Figure  3-2  Insures  that  the  original  passband  shape  from  level  one 
is  maintained  essentially  intact  through  an  arbitrary  number  of  interpolation  levels.  Thus 

I in  the  frequency  domain  the  net  effect  of  each  additional  level  is  essentially  one  of  scaling 

the  passband  portion  of  the  frequency  function  down  by  a factor  of  two. 

The  sampling  rate  equalization  artifice  used  above  in  order  to  make  meaningful  the 
frequency  transfer  function  concept  may  be  avoided  completely.  The  impulse  response 
of  the  cascade  may  be  obtained  by  inputting  a single  unit  value,  followed  by  enough  zero 

I samples  to  force  out  the  entire  impulse  response  sequence.  Once  these  numbers  are  avail- 

able, a normal  and  equivalent  filter  structure  may  be  constructed.  Of  course,  if  this 
equivalent  filter  is  used  for  interpolation  purposes  computational  advantage  should  be  taken 
of  the  fact  that  only  one  in  every  2^  input  samples  is  nonzero.  Tliis  topic  will  be  touched 
upon  again  in  the  paper.  The  immediate  concern,  however,  is  the  class  of  filters  that  may 
be  derived  from  the  impulse  response  of  the  cascade  of  Figure  3-1. 

4.2  DESIGN  EXAMPLES 

Filters  derived  from  the  interpolator  cascade,  not  surprisingly,  have  many  of  the 
same  characteristics  of  the  original  interpolation  process  as  discussed  relative  to 
Figure  2-2  (also  see  Tables  4-1  and  4-2).  Computer  program  C432  in  the  Appendix 
contains  12-selectable  cascade  designs  for  interpolation  Identified  by  parameter  NT. 

Designs  for  stopbands  of  -65,  -50,  -40,  and  -25  dB  are  available.  The  first-level  N value 
selections  for  these  designs  are  shown  circled  and  numbered  in  Figure  2-2.  Each  cascade 


4-3 


M 

§ 

IO 

m 

n 

a* 

IO 

CO 

& 

00 

Cl 

CO 

00 

a* 

S 

II 

a 

•5 

j a* 

00 

a« 

A 

n 

o3 

& 

H 

lO 

$ 

04 

« 

•sf 

olN 

z 2 


u rr- 

& “*.« 


gd- 

A 

•o  Z 

la" 
,,  0 
P-2 

. CJ  (N 

* s 

Q 8 H 

o y 


conhhI^wnh 


N H H 


N H H H 


00  to  T)I  M 

t-  t-  A i-H 

B N H H 


1 

rH 

o 

IO 

Cl 

IO 

a* 

CO 

CO 

35 

S 

19 

3 

A 

to 

00 

M M M 

rH 

00 

IO 

CO 

Cl 

CO 

H 

rH 

CO 

§ 

IO 

H 

H 

tH 

IO 

IO 

IO 

fH 

O 

t- 

IO 

CO 

o> 

l> 

c- 

IO 

CO 

CM 

H 

§3 

H 

13 

c-  1 

0> 

IO 

CO 

u 

z * 

28 

20 

15 

A 

46 

34 

25 

15 

15 

12 

A 

ID 

00 

« a* 

00 

Si 

20 

15 

A 

45 

34 

25 

15 

15 

12 

A 

IO 

•3 

PU  « « 

S 

19 

14 

A 

45 

33 

24 

15 

14 

11 

A 

ID 

Charg 
L = 
L 2 

8 

18 

13 

co 

43 

32 

8 

14 

13 

10 

ao 

IO 

8 

-o* 

o 

to 

ao 

oo 

o 

C| 

o 

oo 

to 

*H 

fH 

n 

A 

Cl 

H 

a _,  tnooioinoommooin 


£ I H^W^lO(Oh®)®OHW 

Z I ^ 1-1  »H 


*si 


TABLE  4-2 

INTERPOLATOR  TYPE  DESCRIPTIONS 


Interpolator 

Type 

(NT) 


1st 
Level 
y (%) 


Stopband 
Level (dB) 


Sampling 

Rate 

Type 

y Progression 
Level  2-5 
(%) 

Critical 

60, 30,15,7.5 

Critical 

60, 30, 15,7.5 

Critical 

60,30, 15,7.5 

Critical 

60.30.15.7.5 

Critical 

60,  30, 15,7.5 

Critical 

60, 30, 15,7.5 

Critical 

60,30, 15,7.5 

Critical 

60.  30.15.7.5 

Oversampled 

30, 15,7.5,3.75 

Oversampled 

30, 15,7.5,3.75 

Oversampled 

30, 15,7.5,3.75 

Overs  ampled 

30, 15,7.5,3.75 

N.  Values 
byKLevel 
1-5 


19, 5, 3,  2, 2 
14,4,  2,  2,1 
10,3,2,1,1 
6.2. 1.1.1 


design  was  obtained  by  drawing  horizontal  lines  on  Figure  2-2  at  the  four  stopband  levels 
(taking  into  account  the  6-dB  gain  provided  by  interleaving).  Safe  N values  at  the  various 
cascade  y lines  (see  Table  4-2  for  the  y-progreslons)  were  selected  for  five  levels  and 
the  appropriate  coefficients  for  these  interpolator  designs  were  incorporated  into  the 
program. 

As  an  example,  the  NT  = 2 processor  provides  for  -50  dB  stopbands  and  starts  with 
a first  level  N of  7 at  a y value  of  80%.  The  second  level  y is  60%  with  an  N value  of  4. 

The  third,  fourth,  and  fifth  levels  have  N values  of  2,  2,  and  1,  and  y values  of  30,  15, 
and  7. 5%.  Designs  NT  = 1-8  ut>e  extended  second-level  y values  and  are  suitable  for  low- 
pass  filter  design.  Designs  NT  = 9-12  are  meant  only  for  Interpolation  in  the  oversampled 
case.  Their  "stopband"  usage  mode  for  the  derived  filters  will  not  be  considered  here. 

If  one  considers  use  of  the  Figure  3-1  cascade  for  the  derivation  of  impulse  response 
data,  it  is  interesting  to  follow  the  development  of  this  response  through  the  cascade.  At 
each  stage  the  input  response  is  stretched  and  new  values  are  interleaved  into  the  vacated 
spaces  as  per  Figure  2-3.  The  passband  region  of  the  frequency  response  function  is  essen- 
tially scaled  down  by  two  to  one  in  the  process.  After  a sufficient  number  of  Interpolator 
stages  are  traversed,  continued  scaling  may  be  done  by  simple  linear  interpolation  as 
previously  discussed.  Hence  the  frequency  scaling  operation  can  be  extended  easily  in 
most  cases  beyond  values  of  practical  interest. 

Examples  of  the  above  techniques  are  demonstrated  in  Figures  4-1  through  4-7.  (All 
frequency  functions  in  this  paper  are  normalized  to  foldover  frequency. ) For  this  demon- 
stration, design  NT  = 5 was  chosen  which  is  admittedly  a "showboat"  design.  NT  = 5 yields 
the  best  performance  (transition  region  width  is  22%  of  passband  width,  W^/f^  = 0.  22,  side- 
lobe  level  = -65  dB),  but  also  has  the  largest  processor  demand.  The  tests  run  were 
very  simple:  the  impulse  responses  from  the  various  levels  of  Figure  3-1  were  individually 
extracted  and  written  on  disc  files  by  the  time-share  program  of  the  Appendix.  At  a later 
time  a batch  program  read  these  files,  performed  on  8192-point  Fast- Fourier  transform  (FFT) 
and  plotted  the  results,  normalized  to  the  peak  level  found  in  the  entire  positive  frequency 
space.  The  plots  show  the  entire  4097-point  spectrum  with  the  foldover  frequency  referenced 
to  unity  at  the  far  right  of  each  plot. 
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Figure  4-3.  Filter  Function,  y = 90%,  NT  = 5,  L = 3,  f = 0. 1125,  Nf  = 171 
Peak  Passband  Ripple  = -61. 95  dB  ** 


Figure  4-5.  Filter  Function,  y.  = 90%,  NT  = 5,  L = 5,  f =0. 028125,  N,  = 693, 
Peak  Passband  Ripple  = -62. 80  dB  **  1 


Figure  4-7.  Filter  Function,  y . = 90%,  NT  = 5,  L = 7,  f =0. 00703125,  N.  = 2775 
Peak  Passband  Ripple  = -61. 26  dB  p 1 


The  impulse  responses  which  result  (see  par.  4. 5)  from  this  synthesis  procedure  and 
computer  program  have  a central  term  of  unity  (the  lone  input  nonzero  sample)  which  is 
delayed  in  the  file  by  the  count  given  by  Equation  (4-1).  The  response  is  symmetric  about 
the  central  term  (linear  phase  delay)  and  the  number  of  active  coefficients  on  each  side  of 
the  central  term,  Np  is  also  given  by  Equation  (4-1). 

Knowledge  of  the  number  of  input  pulses,  N , required  to  fully  load  the  interpolation 

c 

cascade  is  often  useful.  Ibis  number  defines  the  transient  or  "charge-up"  period  which 
precedes  the  normal  steady-state  operation  of  interest.  This  number  may  be  calculated 
from 


Nc  = 


L — 1 


E *<L*1-klNk 


k = 1 


(4-3) 


This  calculation  is  performed  by  the  computer  program  described  in  Section  6,  par.  6.4 
and  tabulated  data  appears  in  Table  4-1. 

The  level  one  filter  frequency  function  is  shown  in  Figure  4-1  and  the  Chebyshev 
behavior  in  both  the  pass  and  stopbands  comes  as  no  surprise.  The  stopband  is  at  about 
the  -67  dB  level  which  checks  very  well  with  the  -61  dB  circled  first-level  choice  for 
NT  = 5 in  Figure  2-2.  The  passband  edge  fp  is  0.45  relative  to  foldover  (one-half  of 
0.90).  The  transition  band  width  is  always  a fixed  22%  of  the  passband  in  these  figures 
due  to  the  frequency  scaling.  Hence  only  the  passband  frequency  value  fp  need  be  given  to 
define  the  filter  behavior.  The  Nf  value  for  this  filter  is  37,  however  18  of  these  coeffi- 
cients are  zero. 
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Figure  4-2  shows  the  level  two  filter  function  with  fp  = 0.  225.  The  level  one  sldelobes 
are  compressed  In  frequency  due  to  the  stretching  in  time.  If  there  were  no  level  two  fill 
of  the  vacated  spaces,  the  level  one  passband  would  appear  at  the  far  right  of  Figure  4-2. 

The  net  effect  of  the  level  two  fill  in  the  frequency  domain  is  very  nearly  a multiplication 
of  the  level  one  passband  with  level  two  stopband.  Thus  the  -70  dB  sldelobes  at  the  right 
of  the  figure  are  the  coarse  sldelobes  of  the  level  two  interpolator  which  has  an  N value  of 
5 (compared  to  19  for  the  first  level).  The  -70  dB  sldelobes  agree  very  well  with  the 
-64  dB,  N = 5,  y =60%  value  from  Figure  2-2. 

Note  also  in  Figure  4-2  that  the  level  one  sldelobes  are  being  affected  by  the  level 
two  transition  region  which  extends  from  0. 3 to  0. 7 and  is  6-dB  down  at  0.  5.  As  additional 
levels  are  traversed,  the  transition  region  effects  cumulate  to  drive  much  of  the  sidelobe 
structure  below  -100  dB  as  will  be  seen.  The  for  this  filter  is  83. 

The  "stretch-and-flll,,  pattern  effects  of  Figure  4-3  for  the  level  three  response 
show  a stopband  at  about  -84  dB.  This  again  agrees  well  with  the  Figure  2-2,  N = 3 value 
at  y - 30%  of  -78  dB.  The  transition  region  for  level  three  extends  from  0. 15  to  0.  85. 

The  passband  frequency  fp  for  this  filter  is  0. 1125  and  the  N^  is  171. 

The  four-level  equivalent  filter  of  Figure  4-4  has  an  f of  0. 05625  and  an  Nj  of  345. 

The  level  four  stopband  at  -82  dB  corresponds  to  a -76. 5 dB  level  at  N = 2 for  y = 15%  in 
Figure  2-2.  The  transition  region  of  level  four  extends  from  0. 075  to  0. 925. 

The  level  five  response  of  Figure  4-5  has  an  f = 0. 028125  and  an  Nj  of  693.  No 
response  is  seen  near  the  foldover  frequency  since  an  N of  2 at  y =7.5%  (level  5 design) 
places  the  stopband  response  near  -107  dB  and  hence  out  of  sight.  The  transition  region 
here  extends  from  0.0375  to  0.9625. 

The  level  six  performance  shown  in  Figure  4-6  is  especially  interesting  as  this 
represents  a departure  in  operating  procedures.  The  interpolation  program  provides 
Figure  2-3  type  processors  at  each  of  the  first  five  levels.  Aftev  the  fifth  level  simple 
linear  Interpolation  is  provided.  The  -59  dB  control  line  (for  -65  dB  performance)  drawn 
in  Figure  2-2  which  was  followed  for  levels  2,  3,  4,  5 at  y values  of  60,  30,  15,  7. 5%  would 
Indicate  one  more  stage  of  N = 1 at  3. 75%  for  the  sixth  level.  The  optimum  N = 1 designs 
have  6-dB  less  peak  error  than  straight  Interpolation  as  was  discussed  relative  to  Equa- 
tions (2-13)  through  (2-16);  linear  interpolation  is  equivalent  to  bj  = 1 in  Equation  (2-17). 

The  N = 1 point  of  Figure  2-1  for  y = 3. 75%  indicates  that  linear  interpolation  at  the  sixth 
level  would  provide  a stopband  level  of  -62  dB.  (The  6-dB  gain  of  the  note  of  Figure  2-1 
is  exactly  offset  by  the  6-dB  loss  of  linear  interpolation. ) 
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The  actual  peak  response  level  noted  near  the  foldover  frequency  of  Figure  4-6  Is 
closer  to  -66  dB  which  Is  within  the  -65  dB  performance  specification.  The  restricted 
response  width  out  of  level  five  (fp  = 0.028125)  is  the  saving  factor.  This  f value,  when  sub- 
stituted into  Equation  (2-13)  for  f/(R/2)  with  ^ = 1,  reveals  an  error  level  of  -60. 21  dB 
at  band  edge  and  when  corrected  by  the  6-dB  gain  from  interleaving,  gives  a -66.  21  dB 
response  level  which  confirms  the  Figure  4-6  results.  Note  also  the  sharpness  of  the 
response  drop  on  the  left  side.  Comparison  of  measurements  from  0 to  passband 
edge  and  from  1 to  the  sharp  drop  in  the  null  region  response  also  confirm  the  causal 
mechanism  cited.  At  this  level  f is  0. 0140625  and  N^  is  1387.  With  linear  interpolation 
or  N = 1 one  might  consider  that  {he  entire  frequency  space  is  a transition  region. 

Figure  4-7  shows  level  seven  output  which  is  the  result  of  five  levels  of  1:2  interpola- 
tion and  one  level  of  1:4  linear  interpolation.  A calculation  similar  to  the  one  for  level  six 
predicts  a stopband  peak  of  -78.  25  dB  which  is  confirmed  by  the  Figure.  The  fp  here  is 

0.  00703125  and  the  Nf  is  2775.  It  is  clear  that  this  scaling  process  could  be  continued 
further  without  difficulty. 

The  three  filter  functions  in  the  figures  which  follow  were  all  taken  at  level  five  in 
the  y = 90%  set  so  that  fp  = 0. 028125  and  W^/fp  = 0.  22  for  all  of  these.  Figure  4-8  shows 
the  NT  = 8 design  which  was  specified  for  -25  dB  stopbands  and  used  an  sequence  of 
6,  2,  1,  1,  1.  This  high  stopband  level  keeps  the  entire  stopband  region  in  view.  The  inter- 
polator levels  responsible  for  sidelobe  groups  seen  in  left-to-rlght  order  are: 

1,  2,  1,  3,  1,  2,  1,  4,  1,  2,  1,  3,  1,  2,  1,  5. 

For  this  high  sidelobe  design  the  mainlobe  of  passband  ripple  is  large  enough  to  be 
visible  in  the  plot.  Calculations  made  with  an  auxiliary  test  program  indicate  that  the  peak 
passband  ripple  deviation  is  at  the  level  of  -22  dB.  This  compares  with  a peak  ripple 
level  of  -26. 7 dB  at  the  output  of  level  one.  The  peak  passband  ripple  level  increased  by 
4. 7 dB  in  transit  through  the  four  levels  two  through  five. 

Figure  4-9  is  next  in  this  series  with  an  NT  value  of  7 which  is  specified  at  -40  dB 
sldelobes.  The  sequence  here  is  10,  3,  2,  1,  1. 

Finally,  for  this  group  Figure  4-10  for  NT  = 6 gives  results  for  the  -50  dB  design 
which  used  an  N^  sequence  of  14,  4,  2,  2,  1.  The  results  are  consistent  with  expectations. 

I I 
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These  plots  of  derived  filter  transfer  functions  have  allowed  the  detailed  workings 
of  this  synthesis  procedure  to  be  examined.  This  data  has  also  confirmed  quantitatively  the 
performance  predicted  by  the  theory  of  the  synthesis  procedure.  Comparisons  of  these 
filters  with  those  of  more  conventional  design  in  terms  of  efficiency  in  application  and  other 
similar  factors  are  given  in  par.  4.4  of  this  section. 

4. 3 DESIGN  PROCEDURES 


The  basic  purpose  of  a filter  is  to  pass  one  band  of  frequencies  and  exclude  another 
band.  One  is  concerned  therefore  with  ripple  in  the  passband,  attenuation  in  the  stopband 
and  the  width  of  the  transition  region.  At  first  level  output,  the  passband  edge  f is 


f = 
P 


2 


(4-4) 


The  stopband  edge  f is 
8 


fs  = 1- 


with  a transition  region  width  wT  of 


(4-5) 


WT  = f8"fp  = 1-?l  <4-6> 

all  relative  to  the  foldover  frequency  which  equals  one-half  the  sampling  rate.  (Sampling 
rate  is  2R  at  the  output  of  level  one. ) 

The  overlap  in  second  level  y ^ required  to  place  transition  regions  in  the  stopband 
requires  that 

> yl 

V2  = 1 “ — (4-7) 

From  this  level  on  the  y values  may  be  halved  at  each  level.  Equation  (4-7)  shows  that 
y j values  of  2/3  or  less  result  In  a y ^ which  is  equal  to  or  greater  than  y ^ which  increases 
filter  length.  At  y j = 2/3  the  transition  region  width  equals  the  passband  bandwidth.  This 
represents  a rather  poor  filter  in  any  event. 

The  Appendix  computer  program  provides  y 1 values  of  80%  and  90%  for  filtering. 

At  80%  the  transition  bandwidth  is  one-half  the  passband  extent  so  that  w^/f  = 50%  for 
y i = 80%.  For  y j = 90%  we  have  wj/fp  * 22%,  a better  filter  but  a costlier  one. 
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The  design  procedure  is  quite  simple: 

1.  Select  transition  region  relative  width  y1  = 80%  (NT  = 1 - 4)  or  y ^ = 90% 
(NT  =5-8),  currently  available. 

2.  Select  stopband  level  (-25,  -40,  -50,  -65  dB),  NT  is  now  Identified.  (See 
Tables  4-2  and  4-2. ) 

3.  Normalize  desired  filter  cutoff  frequency  f by  use  of  Equation  (4-13)  to 

c 

obtain  fp' . Compute  frequency  compression  ratio  Cf 

Cf  = vv  <4 

4.  If  is  acceptably  close  to  a power  of  2,  compute 

L = LEV  = log2  (Cf)  (4 

(if  this  is  not  possible  go  to  step  6. ) 

5.  Use  the  Appendix  program  (XLINT,  LINT)  to  generate  the  desired  filter 
impulse  response.  Task  complete. 

6.  If  the  Cj  factor  cannot  be  approximated  sufficiently  well  by  a power  of  2, 
find  and  L and  J such  that 


P at 

Cf  - 


(4-10) 


(See  Section  5 par.  5.4  for  details  of  L,  J determination. ) 

(a)  Using  LEV  = L get  impulse  response  data  as  in  step  5. 

(b)  Derive  desired  symmetric  filter  Impulse  response  by  using  central  term 
(unity)  of  (a)  above  and  selecting  every  J**1  term  thereafter.  Task  complete. 

For  an  example  of  a "J-derived"  filter  consider  the  following:  Let  sharpness  of  cutoff 

require  y^  = 90%.  Let  sldelobe  level  specification  be  < -60  dB,  so  NT  = 5 is  selected. 

Assume  clock  period  in  processor  hardware  of  200  ps  or  5000/s  repetition  rate,  and  assume 

filter  cutoff  frequency  f = 225  Hz  ± 2%.  Normalize  this  cutoff  frequency  to  the  foldover 

c 

frequency  as  per  Equation  (4-13). 


V - 


= 0.09 
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then  from  Equation  (4-6) 


Ct  = y./f’  = 0.9/0.09  = 10 
I ' 1 p 

This  is  not  a power  of  two,  go  to  step  6 and  try 
~ 27 

Cf  = 10  = 13  = 9‘8462 

This  will  give  f = 228. 52  Hz  which  is  within  the  ±2%  tolerance, 
c 

Use  LEV  = 7,  NT  = 5 and  get  impulse  response  from  Appendix  program  C432.  From 
Equation  (4-1)  we  find  that  the  central  peak  is  delayed  by 

Delay  = 2775  + 1 = 2776. 

(The  central  peak  is  in  cell  number  2777  of  the  array. ) Equation  (4-1)  states  that  there  are 
Nf  = 2775 

active  coefficients  on  each  side  of  unity.  After  every  13th  element  is  selected  and  the  new 
impulse  response  is  created,  it  will  have  a delay  of  213  and  anNj  = 213.  Figure  4-11  shows 
the  frequency  function  for  this  filter.  The  impulse  response  involved  in  Figure  4-7  was 
censored  to  yield  the  impulse  response  of  the  Figure  4-11  data. 

Figure  4-11  could  also  be  compared  to  Figures  4-3  and  4-4  which  are  the  closest 
"normal-sequence"  filters.  Apparently  no  real  violence  was  done  to  the  basic  frequency 
response  by  the  resampling  process.  This  suggests  that  very-high  LEV  Impulse  responses 
be  calculated  and  stored  for  a particular  combination  of  transition  width/passband  ratios 
and  sidelobe  levels.  These  stored  responses  could  then  simply  be  resampled  by  some 
appropriate  factor  J for  specific  cutoff  frequencies. 

In  this  regard  the  three  equations  which  follow  may  be  helpful  in  defining  the  filter. 
Once  the  transition  bandwidth  to  passband  frequency  ratio  (w^/f^)  is  known,  the  first- 
level  y becomes 


(4-11) 
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Figure  4-11.  Filter  Function  for  "J -Derived"  Filter 


A specified  cutoff  frequency  f is  normalized  to  its  foldover  frequency  to  yield  f ' 

c p 


f ' = f /(one-half  sampling  rate)  (4-13) 

P ® 

which  may  then  be  substituted  into  Equation  (4-8). 

4.4  COMPARISON  TO  CONVENTIONAL  DESIGNS 

A rule-of-thumb  guide  often  used  to  estimate  the  total  length  K of  a low-pass  digital 
filter  is  (Ref.  4)  given  by  Kaiser  as 


K = 


-!0  log1Q  (El  - 15 
14  (W) 


(4-14) 


where  Et,  E0  are  pass  and  stopband  peak  ripple  values  and  W is  the  transition  region 
width  expressed  relative  to  sampling  rate.  When  expressed  in  terms  of  filter  half  length 
Nf  and  transition  width  wT  relative  to  the  foldover  frequency  one  obtains 


Nf  " 


-10  log10  (Ej  Eg)  - 15 
14wl 


(4-15) 


The  filter  synthesis  procedure  described  in  this  paper  provides  no  direct  control 
of  passband  ripple.  The  individual  level  processors  provide  a passband  gain  proportional 
to  (1  + 6 ) where  6 is  an  error  directly  related  to  stopband  level.  Thus,  through  two 
stages,  the  passband  gain  is 


(1+5.)  (1  + 5J  = 1 + «.  +6, 


(4-16) 
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All  stages  have  nearly  the  same  6 value  and  since  the  stages  have  different  y values,  it 
seems  reasonable  to  assume  that  the  <5  's  will  add  on  an  rms  basis.  A fair  estimate  of 
ripple  might  then  be  obtained  and  substituted  into  Equation  (4-15)  so  that  equivalent 
standard-design  filter  sizes  might  be  estimated.  This  yields 


„ -10log10(^L  62)  -15 

Nf  = 

* -20  log1Q  (<5 ) - 5 log1Q  (L)  - 15 
Nf  = 


(4-17) 


(4-18) 


The  first  numerator  term  above  is  the  negative  of  the  stopband  level  in  dB. 


The  filters  described  in  this  paper  may  be  sized  assuming  that  w^/fp  is  given  so  that 
may  be  computed  from  Equation  (4-11).  Convenience  is  gained  with  no  loss  of  generality 
if  comparisons  are  normalized  to  level  one  output  conditions.  This  results  in  a w^  value  of 
(1  - Yj)  according  to  Equation  (4-6).  The  wT  at  level  L,  wTL  would  then  be 


(i  - yx) 

WTL  " 2L  - 1 


(4-19) 


This  allows  Equation  (4-18)  to  be  written  as 


-20  log1Q  (<5 ) - 15  - 5 log1Q  (L) 
14  (1  - yy) 


- 1 


(4-20) 


Once  the  L value  is  selected,  Equation  (4-1)  or  Table  4-1  may  be  used  to  obtain  the 
required  Nj  which  may  be  compared  to  that  obtained  from  Equation  (4-20). 

Letting  pN^  be  the  coefficient  count  of  the  filter  design  of  this  paper  and  be  the 
corresponding  Kaiser  (Ref.  4)  estimate,  the  results  of  Figure  4-12  obtain.  The  paper 
design  appears  to  require  more  coefficients,  but  the  ratios  are  reasonable.  The  passband 
ripple  problem  of  the  paper  procedure  shows  to  its  disadvantage  at  the  lower  stopband 
suppression  levels,  as  expected.  For  the  higher-quality  filters  the  paper  design  compares 
reasonable  well  with  more  conventional  synthesis  procedures.  (Later  passband  ripple 
studies  indicate  that  the  ratios  of  Figure  4-12  are  too  large  for  the  higher  L values). 
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4.  5 TIME  DOMAIN  CONSIDERATIONS 


In  previous  discussions  the  filter  frequency  function  was  examined  at  different  levels 
in  Figure  3-1  structure.  A similar  study  will  be  made  here  of  the  filter  impulse  response 
function.  The  case  NT  = 5 interpolator  places  serious  restrictions  on  stopband  level  and  the 
processor  is  nontrivial  at  all  levels  ( Nk  = 19,  5, 3,  2, 2).  Figure  4-13  shows  the  first  level 
impulse  response  which  has  an  N^  = 37.  (The  sample  values  are  the  significant  items 
of  concern,  theCalComp  plotter  draws  straight  lines  between  points.  The  discontinuities 
in  slope  indicate  point  location. ) The  frequency  function  for  this  data  appears  in 
Figure  4-1. 

The  second  level  processor  is  also  a sophisticated  one  (N2  = 5)  so  that  when  the  first 
level  input  values  are  interleaved  with  new  derived  values  the  output  of  Figure  4-14  results. 
Considerable  smoothing  has  been  done  but  slope  breaks  are  still  evident.  The  N^  value 
is  83.  The  corresponding  frequency  function  is  that  of  Figure  4-2. 

The  third  level  processor  has  an  Ng  = 3 and  its  "stretch-and-flll"  action  produces  the 
considerably  smoothed  data  of  Figure  4-15  (frequency  function  is  Figure  4-3).  Slope 
granularity  is  still  in  evidence  near  local  peaks.  The  for  this  data  is  171. 

Fourth  level  output  in  Figure  4-16  is  very  smooth;  the  processor  sophistication  has 
dropped  to  N^  = 2 here.  N^  is  345  and  the  corresponding  frequency  function  is  that  of 
Figure  4-4. 

Another  Ng  = 2 stage  at  level  five  produces  the  output  shown  in  Figure  4 17  which 
has  an  Nj  = 693  and  corresponds  to  the  frequency  function  of  Figure  4-5. 

Beyond  this  point  in  the  chain,  linear  interpolation  may  be  used.  The  "J-derived" 
filter  of  Figure  4-11  has  the  impulse  response  shown  in  Figure  4-18  where  N^  = 213.  This 
filter  was  derived  by  using  L = 7 and  then  selecting  every  13th  sample  of  the  resulting 
impulse  response. 

This  rather  attractive  smoothing  of  the  impulse  response  with  an  increase  in  L does 
not  always  occur.  Consider  Figures  4-19  and  4-20.  Here  NT  = 8 (y  1 = 90%,  -25  dB 
stopband)  and  the  {N^}  sequence  is  6,  2,  1,  1,  1.  Because  of  the  high  sidelobe  levels,  the 
stages  of  Figure  3-1  beyond  the  second  are  just  modified  linear  Interpolators  (N^  = 1). 

The  level  one  output  is  coarse  as  expected  (N^  = 11).  However  note  that  the  L = 7 output 
of  Figure  4-20  shows  considerable  slope  granularity.  The  high  sidelobe  specification 
(-25  dB)  permits  this  kind  of  crude  "stretch-and-flll"  operation.  (The  Y-axis  scaling  for 
all  impulse  response  plots  is  the  same.  There  are  subtle  differences  in  the  X-axis 
positioning  and  scaling  for  the  different  plots. ) 
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Figure  4-14.  Filter  Impulse  Response,  NT  = 5, 
Frequency  Function  is  Figure  4-2 


Figure  4-15.  Filter  Impulse  Response,  NT  = 5,  LEV  =3,  = 171. 

Frequency  Function  is  Figure  4-3 


Figure  4-17.  Filter  Impulse  Response,  NT  = 5,  LEV  =5,  N,  = 693. 
Frequency  Function  is  Figure  4-5  1 


Figure  4-19.  Filter  Impulse  Response, 


4.6  PASSBAND  RIPPLE  CONSIDERATIONS 


Small  deviations  from  unity  gain  in  the  passband  at  each  level  result  in  a deviation 
from  unity  in-band  gain  of  the  resultant  derived  filter.  The  derived  filter  passband  gain 
deviation  is  approximately  equal  to  the  sum  of  the  individual  level  deviations.  (The  derived 
filter  gain  must  be  normalized  as  the  nominal  voltage  gain  is  2L  in  the  passband. ) In  this 
design  procedure  the  stopband  peak  levels  are  under  strict  control.  The  passband  peak 
ripple  levels  are  not  similarly  constrained.  The  nature  of  passband  gain  error  cumulation 
is  somewhat  complex.  Measurements  of  peak  in-band  ripple  level  for  some  of  the  derived 
filters  discussed  here  have  been  made  with  rather  encouraging  results. 

The  worst-case  assumption  of  coherent  addition  of  the  deviations  was  found  to  be 
far  too  pessimistic.  The  rms-addition  assumption  of  Section  4,  par.  4.4  was  found  to 
slightly  understate  the  peak  ripple  level  in  the  L = 2,  3 region.  For  higher  L values  the 
rms  assumption  becomes  too  pessimistic.  For  example  in  NT  = 5,  L = 7 case  the  peak 
in-band  ripple  is  about  5. 8 dB  worse  than  at  the  output  of  level  one.  An  rms  addition  assump- 
tion would  predict  an  8. 5 dB  degradation. 

If  in-band  ripple  is  of  very  serious  concern  one  can  p roceed  with  the  designs  as  out- 
lined previously.  When  complete,  the  precise  normalized  gain  of  the  derived  filter  may  be 
calculated.  If  this  gain  is  represented  by  1 + 6 (f)  then  the  compensating  gain  function 
1 - 6 (f)  may  be  immediately  specified.  The  procedures  of  Section  2,  par.  2.  2 may  then 
be  used  to  determine  the  coefficients  of  a pre-emphasis  correction  filter.  For  the  specific 
case  of  NT  = 5,  L = 7 the  writer  used  this  procedure  to  reduce  peak  in-band  ripple  to  level 
one  (stopband)  values.  The  pre-emphasis  filter  complexity  required  was  equal  to  that  of 
the  first  level  filter  the  original  cascade  (N^  = 19).  An  alternative  approach  might  be 
to  overspecify  stopband  level  at  the  outset  to  provide  room  for  in-band  ripple  growth 
through  the  cascade. 

In  order  to  keep  the  ripple  problem  in  perspective  Figure  4-21  shows  the  variation  of 
pass  band  gain  with  ripple  magnitude.  Ripple  levels  in  the  -20  dB  region  cause  less  than 
a 1-dB  gain  variation.  Gain  variations  are  the  order  of  0. 05  dB  at  -40  dB  ripple  levels. 

In  many  engineering  applications  the  sidelobe  specification  will  result  in  ripple  levels 
which  produce  negligible  pass  band  gain  variations  in  spite  of  ripple  growth  through 
the  cascade. 
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PASSBAND  GAIN  VARIATION  (dB) 


Figure  4-21.  Baseband  Gain  Variation  as  a Function  of  Ripple  Value 
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SECTION  V 
APPLICATIONS 


5.1  MODE  CONVERSION.  COMPLEX-TO-REAL 

Those  who  do  practical  digital  signal  processing  work  usually  deal  with  a variety  of 


processing  software  and  hardware  modules.  One  must  also  deal  with  a number  of  different 
formats  in  the  data  supplied  as  input  or  specified  for  the  processed  output  data.  One 
of  the  most  significant  format  characteristics  concerns  the  nature  of  the  sampling  process 
Itself.  The  individual  sampling  operation  may  involve  one  real  number  (a  "real"  sample) 
or  an  ordered  pair  of  real  numbers  (a  "complex"  sample). 


Consider  the  low-pass  signal  spectrum  illustrated  in  Figure  5-l(a).  The  signal  shown 


has  an  upper  frequency  limit  of  2f  so  that  the  sampling  rate  R must  be  4 f if  the  usual 

c c 


single-number-per-sample  ("real")  low-pass  sampling  is  done.  It  often  happens,  however, 
that  the  signal  band  of  interest  is  centered  at  a frequency  which  is  high  compared  to  the 
bandwidth  of  the  signal.  In  this  case  sampling  at  a rate  equal  to  twice  the  highest  frequency 
is  very  inefficient  and  one  of  two  alternate  procedures  is  normally  followed. 


For  purposes  of  illustration  assume  that  the  bandpass  signal  has  a bandwidth  of  2f  . 

c 


The  signal  could  be  bandpass  filtered  then  frequency-shifted  down  to  a new  center  frequency, 


f , as  shown  in  Figure  5-l(b).  At  this  point  the  usual  low-pass  (real)  sampling  could  be 
c 


made  at  rate  4 f . This  procedure  requires  a front-end  filter  tailored  to  the  bandpass 

C 


signal  spectrum,  a mixer,  local  oscillator,  and  post-mixer  low-pass  filter  as  shown 
in  Figure  5-l(b). 

An  alternate  approach,  shown  in  Figure  5-  1(c),  uses  product  detection  by  cosine  and 
sine  signals  from  a local  oscillator  at  the  band  center  frequency,  followed  by  low-pass  filters 
to  produce  the  familiar  i (t)  and  q (t)  signals.  (In  the  general  case  these  are  not  Hilbert 
transform  related. ) These  in-phase  and  quadrature  signal  components  have  a high  frequency 


of  fc>  These  signals  are  then  each  sampled  at  a rate  2ffi  to  produce,  at  each  sample  time, 


an  ordered  pair  of  real  numbers  (i^,  qk ).  This  pair  is  often  referred  to  as  a "complex" 
sample. 


These  two  sampling  techniques  are  theoretically  equivalent.  Each  has  its  own 
practical  problems  and  advantages  and  both  are  commonly  encountered  in  practice.  Data 
processing  modules  are  not  always  designed  to  accept  or  produce  either  format.  Hence 
one  aften  encounters  the  practical  need  for  mode  conversion  from  real  to  complex  or  from 
complex  to  real.  The  complex-to-real  conversion  will  be  discussed  first  with  the  aid 
of  Figure  5-1. 
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b.  Frequency  Downshift  Followed  by  Baseband  Sampling 


c.  Complex  (Bandpass  Sampling) 


Figure  5-1.  Sampling  Techniques 


It  may  be  seen  from  Figure  5-l(c)  that  three  quantities  define  the  bandpass  signal: 

the  center  frequency,  i(t),  and  q(t).  If  frequency  translation  of  the  signal  is  properly 

done,  i (t)  and  q (t)  remain  unchanged  and  only  the  center  frequency  changes.  Because  of 

this  Invariance  it  may  be  assumed  that  the  bandpass  signal  has  been  shifted  to  the  center 

frequency  f of  Figure  5-l(a).  The  resulting  time  function  may  then  be  expressed  both  in 
c 

terms  of  real  samples  and  complex  samples  involving  1 (t)  and  q (t)  (which  are  invariant 
with  frequency  translation).  Relationships  between  real  and  complex  sampling  may  therefore 
be  established  and  used  in  mode-conversion  operations. 

The  signal  function  s (t)  may  be  expressed  as 

s (t)  = i (t)  cos  2 ir  ft  - q (t)  sin  2ir  ft  (5-1) 

c c 

A sampling  rate  of  R = 4 fc,  will  yield  time  values  of  t^  = k/4  fc  so  that  Equation  (5-1) 
becomes 


sk  = ifc  cos  — k - qk  sin  y k . 
which  means  that 


(5-2) 


sk  = (-D  ik. 


k = 0,  2,  4,  6, 


(k-1) 

8k  = (~  1)  Qjj  i k = 1,  3,  5,  7, 


(5-3) 


It  may  be  seen  that  the  real  samples  sk  alternately  reflect  the  i (t)  and  q (t)  component 
values.  Complex  samples  are  taken  at  half  the  rate,  but  i (t)  and  q(t)  are  sampled 
simultaneously.  If  pk  represents  a complex  sample  as  per  Figure  5-l(c)  then  it 
follows  that 


(«k>  - - Qj.  " *2.  Qj.  l4.  - 15.  - V V 


w - 


1Q»  » *2*  ” * *4’ 


(5-4) 


ig.  ~ . 

V - » <12*  “»V  " V " ’ 

If  complex  samples  (Pk)  are  given  and  real  sk  samples  are  desired,  Equation  (5-4) 
shows  that  the  even  samples  of  sk  are  directly  available  from  the  ik  terms  of  pk>  The 
odd  samples  of  sk  are  samples  of  q (t)  but  at  times  intermediate  to  the  sampling  times  of 
the  pk  number  pair. 
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A single  level  of  interpolation  of  the  q terms  of  the  complex  samples  is  required. 

Only  the  interpolated  values  are  to  be  kept,  the  original  input  q terms  are  discarded.  The 
delay  inherent  in  the  interpolation  operation  on  the  q sequence  must  be  applied  to  the 
terms,  and  the  sign  changes  in  the  s^  series  of  Equation  (5-4)  must  also  be  accounted  for. 

A subroutine  for  complex-to-real  conversion  based  on  the  above  principles  appears  in  the 
Appendix  as  part  of  C432  and  is  also  described  in  Section  VI. 

5.2  MODE  CONVERSION.  REAL  TO  COMPLEX 

In  Equation  (5-4)  it  is  now  assumed  that  the  sequence  1 s^  } is  given  and  that  the 
sequence  {p^}  is  desired.  For  the  sake  of  variety  only,  synthesis  of  the  complex  sequence 
at  odd-numbered  times  will  be  considered.  The  q portions  of  the  complex  pair  are  available 
directly  from  the  odd  sk  terms  as  shown  in  Equation  (5-4).  The  odd-numbered  i terms 
must  be  generated  by  single-level  interpolation  from  the  even-numbered  i samples.  Only 
the  interpolation  data  is  kept,  the  input  i data  is  discarded.  Delay  compensation  of  the  q 
data  to  match  the  interpolator  delay  must  be  provided,  a1  ig  vith  proper  sign-change 
accounting.  A subroutine  for  real-to-complex  conversion  ' on  the  above  principles 
appears  in  the  Appendix  and  is  described  in  Section  VI.  Fc  ‘ rent  approach  to  the 
real-to-complex  conversion  see  reference  11. 

Conversions  of  the  types  discussed  above  involve  an  Inherent  phase  ambiguity  which 
is  an  integral  multiple  of  n/4. 

5.3  COMPARISON  OF  INTERPOLATION  METHODS 

Schafer  and  Rabiner  (reference  3)  discuss  in  detail  interpolation  using  digital  low-pass 
or  stopband  filters.  In  Section  C of  the  above  reference,  comparisons  to  classical  interpola- 
tion methods  such  as  those  of  LaGrange  are  examined.  These  writers  conclude  that  the 
frequency  response  characteristics  of  the  classical  interpolators  leave  much  to  be  desired. 
They  further  conclude  that  for  most  digital  signal  processing  work  low-pass  or  bandstop 
filters  are  to  be  preferred  over  classical  interpolation  algorithms. 

Attention  can  then  turn  to  a comparison  of  the  iterative  technique  of  Figure  3-1  and 
the  use  of  low-pass  filters  for  interpolation.  It  has  been  shown  in  Section  IV  that  the  filters 
derived  from  the  interpolation  scheme  of  Figure  3-1  appears  to  be  reasonably  competitive 
with  conventional  filter  designs  in  terms  of  processor  demand.  It  is  then  possible  to 
estimate  the  relative  performance  of  the  Iterative  interpolation  method  and  the  filtering 
method  by  use  of  Equations  (4-1)  and  (4-2). 


k 


5-4 


The  O.  result  of  Equation  (4-2)  gives  the  number  of  operations  per  input  value  which 

I* 

produces  2 output  values.  The  Nf  result  of  Equation  (4-1)  gives  the  number  of  filter 

* L 

operations  per  output  value.  However,  In  interpolation  service,  only  one  in  every  2 input 
samples  is  nonzero  so  that  the  computational  load  is  reduced  by  the  factor  2L.  If  the 
comparison  is  based  on  operations  per  output  value,  both  Oj  and  would  be  divided  by 
2L  and  their  ratio  would  remain  unchanged.  Therefore  the  ratio 

L 

2 (L  + 1 - k)  _ (2l  _ i ) 

(5-5) 


provides  the  comparison  desired. 

Figure  5-2  shows  the  plotted  results  of  Equation  (5-5)  for  the  12  Interpolator  designs 
of  the  program  in  the  Appendix.  It  appears  that  the  minimum  processor  demand  advantage 
of  the  iterative  interpolator  is  approximately  a factor  of  two.  As  the  sampling  rate  multi- 
plication factor  increases  the  Iterative  interpolator  advantage  also  Increases.  For  the 
higher-performance  interpolators,  this  advantage  Increase  is  significant. 

(It  could  be  argued  that  for  filters  derived  specifically  by  the  procedures  of  this  paper 
Equation  (5-5)  should  be  multiplied  by  the  factor  (2L  - 1)/2L.  This  results  from  the  exact 
nulls  in  the  Impulse  response  which  have  a spacing  2L.  This  would  bring  the  two  procedures 
into  near  parity  for  low  L values.  The  high  L value  results  of  Equation  (5-5)  however  would 
remain  unchanged.  The  degree  to  which  other  filter  designs  could  take  practical  advantage 
of  impulse-response  nulls  is  difficult  to  assess  in  general.  A heuristic  explanation  of  the 
advantage  of  the  Figure  3-1  iterative  Interpolator  concerns  variations  of  processor  com- 
plexity and  data  rate  through  the  cascade.  Early  in  the  cascade  the  processors  are  the  most 
extensive,  but  the  data  rate  is  lowest.  Later  in  the  cascade  when  the  data  rate  is  high, 
the  processors  are  short. ) The  advantages  of  the  cascade  over  the  lowpass  filter  for 
Interpolation  service  are  also  discussed  in  references  8 and  10. 
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L - LOGj  (SAMPLING  RATE  MULTIPLICATION  FACTOR) 
S - STOP  BAND  LEVEL  (dB) 

N(  - REQUIRED  FILTER  OPERATIONS 

Oj  - REQUIRED  INTERPOLATOR  OPERATIONS 


Figure  5-2.  Processor  Demand  Comparisons,  Lowpass  Filter  vs 
Iterative  Interpolator 


Table  4-1  gives  some  parameters  of  the  interpolators  and  derived  filters  for  levels 
from  1-5.  The  delay  and  tabulations  may  be  used  as  an  aid  in  locating  coefficients  in 
memory  when  the  Appendix  program  is  used  to  provide  equivalent  filter  impulse  response. 

5.4  ARBITRARY  SAMPLING  RATE  MULTIPLICATION  FACTORS 

It  often  occurs,  of  course,  that  the  sampling  rate  multiplication  factor  K is  other 
than  a power  of  two.  The  Appendix  program,  as  written,  will  only  give  powers  of  two 
for  K.  There  are  at  least  two  solutions  to  this  problem. 

The  first  approach  is  that  used  in  the  "J-derived"  filter  of  Section  IV.  The  program 
is  run  with  some  selected  level  L to  yield  a multiplication  factor  of  2L.  The  resulting  output 
data  Is  then  resampled  with  period  J.  The  problem  becomes  one  of  selecting  L,  J so  that 

K = 2L/J  (5-6) 

A simple  solution  involves  selecting  a trial  L and  solving  for  J according  to 


unless 


r 2l 
- n -L-fr-J 


L j?  N + 11  _2K>0 

1 N (N  + 1)  * K 


J = N + 1 


where  [ ] represents  the  Integer  part  of  the  enclosed  quantity.  This  trial  L,  J will  yield 

A 

an  approximate  multiplication  factor9  K,  of 


A 9L 

K --r-  <6- 

An  error  measure  A^  may  be  defined 

A 1 

. A K - K 2^ 

ak  ■ — — nr  - 1 <5 

One  starts  with  the  smallest  L value  which  is  at  least  as  large  as  logg  (K).  Equations 
(5-7),  (5-8)  and  (5-9)  are  used  and  the  resulting  AK  Is  noted.  If  the  error  is  too  large, 
Increment  L by  one  and  repeat  until  the  error  reaches  an  acceptable  level.  An  HP-67 
program  for  Implementing  this  procedure  is  given  in  the  Appendix , and  labeled 
"PROGRAM  1-11-79". 


For  example  let  K = 13.  The  above  procedure  yields  a solution  L = 9,  J = 39  for  1% 

9 10 

error  (2  /39  = 13. 13),  or  L = 10,  J = 79  for  0. 3%  error  (2  /79  = 12. 96).  This  approach 

will  always  yield  the  minimum  L required  to  meet  any  output  sampling  rate  multiplication 

factor  error  specification. 

This  technique  will  work  very  nicely  also  when  the  desired  output  sampling  rate  is 
not  an  integer  multiple  of  the  input  sampling  rate.  For  example  let  a sampling  rate  multipli- 
cation factor  K = ir  be  required  within  0.02%.  The  above  procedure  quickly  yields  L = 9, 

J = 163  (29/163  = 3. 1411,  ^ = 0. 0155%). 

The  second  solution  to  this  problem  involves  a modification  of  the  computer  program. 
One  must  traverse  enough  processors  of  the  Figure  2-3  type  until  the  sampling  rate  is 
sufficiently  high  so  as  to  permit  linear  interpolation.  Once  this  stage  is  reached  the 
sampling  rate  multiplication  factor  at  the  linear  interpolation  level  can  be  any  integer  value. 
Due  to  a foolish  consistency  (hopefully  not  that  of  which  Emerson  spoke)  the  writer  continued 
the  power  of  two  logic  through  the  last  "stage"  when  coding  the  Appendix  program. 

Hindsight  now  makes  it  clear  that  this  stage  should  have  been  coded  for  any  desired 
multiplication  factor,  say  J.  The  overall  sampling  rate  Increase  could  then  be  2^x  J 
where  L is  the  number  of  Figure  2-3  doubler  stages.  The  output  data  could  then  be  censored 
by  a factor  of  2L  to  give  precisely  a 1:J  overall  factor.  This  feature  would  also  lend  more 
flexibility  to  the  operations  discussed  relative  to  Equations  (5-6)  through  (5-9). 

5.5  SAMPLING  RATE  DIVISION 

The  advantages  of  the  cascade  interpolation  technique  of  Figure  3-1  over  the  classical 
low-pass  filter  processing  of  zero-filled  input  data  have  been  discussed  in  detail.  If  one 
studies  Figure  3-1  and  visualizes  2^  output  samples  for  each  input  sample,  simple  faith  in 
the  consistency  of  nature  raises,  on  first  impulse,  the  possibility  of  putting  2^  samples  in 
at  the  bottom  and  getting  one  sample  out  at  the  top.  Is  there  a cascade  equivalent  of  the 
conventional  low-pass  filter/data  censor  combination?  The  answer  is  yes  and,  furthermore, 
the  significant  advantages  of  the  cascade  over  the  low-pass  filter  method  noted  for  inter- 
polation are  again  to  be  enjoyed  in  the  sampling  rate  division  operation,  a FORTRAN  pro- 
gram for  data  rate  division  also  appears  in  the  Appendix  under  deck  name  C434. 

(The  advantages  of  multistage  data  rate  reduction  (decimation)  are  discussed  in 
references  6,  8 and  9.  In  the  case  where  lowpass  filtering  is  to  be  done  while  maintaining 
the  sampling  rate,  reference  9 demonstrates  significant  gains  for  a joint  use  of  rate  reduc- 
tion followed  by  rate  multiplication. ) 
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The  sampling  rate  division  cascade  uses  the  same  level  processor-equivalent 
filters  as  the  interpolation  cascade  but  with  some  significant  differences.  The  processors 
are  used  in  reverse  order.  The  exit  is  always  from  the  level  one  processor,  and  the 
entrance  depends  on  the  number  of  levels  or  L value  of  the  rate-reduction  operation.  For 
example,  If  a 16:1  rate  reduction  Is  specified,  L = 4 and  the  cascade  would  be  entered  at 
the  level  four  design  where  a 2:1  rate  reduction  would  be  achieved  after  filtering.  This 
output  would  then  go  to  a level  3 design  followed  by  a level  2 design,  and  finally  to  the  level 
1 design  for  pre-output  processing.  If  L exceeds  5,  a linear  interpolation  processor  leads 
off  the  chain,  but  this  is  a special  case  which  will  be  treated  separately  later. 

Figure  5-3  illustrates  the  process  at  each  normal  level.  In  Figure  5-3(a)  the  input 
data  at  rate  R is  passed  through  a low-pass  filter  whose  output  is  resampled  for  a 2:1  rate 
reduction.  As  discussed  earlier,  these  filters  have  a unit  central  weight  with  symmetrical 
filter  coefficients  and  interleaved  zero  weights  as  shown.  If  there  are  N coefficient  values 
for  this  filter,  (4N-1)  storage  locations  are  required.  An  exact  equivalent  filter/censor 
combination  requiring  only  (3N+1)  locations  is  shown  in  Figure  5-3(b).  In  this  configura- 
tion the  input  values  are  switched  between  what  turns  out  to  be  the  interpolation  processor 
and  a delay  element.  The  censored  output  is  available  directly  from  the  adder  as  shown. 

It  is  this  latter  processor  that  is  coded  in  0434. 

If  L exceeds  5 the  lead-off  processor  is  a linear  interpolator  followed  by  censoring. 
The  impulse  response  of  the  linear  interpolator  has  a unit  central  term  flanked  by  descend- 
ing sized  terms  equal  to  M/2L~5  where  M = 1,  2,  — (2L~5  -1).  These  terms  are  con- 
tiguous. Following  the  interpolator,  every  2L~5th  value  is  selected  for  output.  This  is 
precisely  equivalent  to  a chain  of  L-5  normal  cascade  processors  having  = 1 and 
ax  = 1/2. 

Sampling  rate  division  at  each  level  is  achieved  by  the  equivalent  of  low-pass  filtering 
followed  by  2:1  data  rate  reduction.  One  final  output  value  is  generated  for  each  2L  input 
values.  The  sampling  rate  is  highest  where  the  processors  are  shortest.  The  processor 
lengths  are  greatest  where  the  sampling  rate  is  lowest  (near  the  output).  The  number 
of  operations  per  output  value  for  the  cascade  divider  is 


(5-10) 


which  is  precisely  the  same  as  the  interpolation  result  O.  of  Equation  (4-2)  which  was 

T * 

based  on  one  input  pulse  yielding  2 output  values. 
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The  delay  through  the  cascade  divider  measured  In  Input  pulses  Is  (see  also 


Equation  (4-1)). 

L 

E 


k=l 


Nk  zL  k+1 


(5-11) 


whether  dealing  with  Interpolation  or  rate  reduction,  the  equivalent  overall  low-pass  filter 
is  the  same  for  a given  (NT,  LEV)  pair.  It  Is  not  surprising  then  to  find  that  the  ratio  of 
operations  required  for  filtering  divided  by  the  number  of  cascade  operations  per  output 
data  value  is,  for  rate  reduction, 


nL-k+1 


£ Nfc  2"  - (2"  - 1) 


k=l 


E Nt2 


k - 1 


k=l 
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which  Is  exactly  the  result  of  Equation  (5-5)  for  interpolation.  It  follows  therefore  that  the 
considerable-computational  efficienty  advantages  of  the  cascade  over  the  filter  shown  in 
Figure  5-2  for  interpolation  apply  here  also.  Note  also  that  the  only  limit  to  rate  reduction 
factor  in  the  present  program  is  the  2^  memory  size  required  of  the  input  array. 

The  impulse  response  determination  of  the  equivalent  low-pass  filter  obtained  by 
driving  the  interpolation  program  C432  with  a unit  value  followed  by  a string  of  zeros  will 
not  work  with  the  rate  reduction  program  C434.  What  one  obtains  in  the  latter  case  is  a 
2L-censored  version  of  the  relevant  impulse  response;  not  a very  useful  result.  This  is 
perhaps  a manifestation  of  the  reversibility  of  the  interpolation  process  as  compared  to  the 
irreversibility  of  the  rate -reduction  process.  The  latter  discards  information,  the  former 
does  not. 

Cascade-derived  impulse  responses  were  used  successfully  and  correctly  in  the 
interpolation  case,  but  care  must  always  be  exercised  because  of  the  sampling  rate  changes 
through  the  cascade.  Even  in  C432  the  Impulse  response  of  a two-level  cascade  is  not  the 
normal  convolution  of  the  individual  level  impulse  responses.  If  one  is  intersted  in  equiv- 
alent filter  Impulse  responses  program  C432  must  be  used;  C434  cannot  be  employed  for 
such  a purpose. 
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The  low-pass  filter  equivalents  of  the  individual  level  processors  have  very  useful 
spectral  properties  for  the  two-to-one  rate  reduction  achieved  in  each  stage.  In  our  frame 
of  reference,  frequency  extends  from  zero  to  unity  (referenced  to  R/2  or  half  the  sampling 
rate).  The  symmetric  coincidence  of  passband  and  stopband  about  f = 1/2  for  these  filters 
is  fortuitous.  The  spectral  components  in  the  stopband  region  will  be  folded  into  the  pass- 
band.  Hence  the  region  of  high  attenuation  is  precisely  where  it  will  do  the  most  good. 

The  spectral  components  near  frequency  f = 1/2  are  attenuated  by  about  6 dB  (trans- 
ition region  center)  and  are  moved  up  to  f = 1 for  the  next  stage.  In  the  next  stage  stopband 
attenuation  is  applied  to  this  reduced  level  aliasing  band  before  being  folded  into  the  pass- 
band.  The  cumulative  attenuation  effects  of  the  transition  regions  as  the  cascade  is 
traversed  keep  the  passband  noise  level  growth  well  under  control  for  high  L values. 


Design  of  a rate  division  processor  is  quite  simple.  The  equivalent  filter  transition 
region  width  and  stopband  attenuation  are  defined  by  the  NT  choice.  The  rate-reduction 
factor  2L  (via  specification  of  LEV  in  C434)  may  be  thought  of  as  defining  a filter  at  input 
level  having  passband  edge  f = y1/2L  where  y ^ refers  to  the  interpolator  level  one  proces- 
sor (which  is  now  at  the  end  of  the  rate-division  chain).  At  the  final  cascade  output, 
censoring  will  have  moved  f to  the  value  y y 

As  a test  of  the  C434  concept,  programs  LINT  and  LDIV  were  run  back-to-back  for 
NT  = 5 and  LEV  = 7.  That  is,  a 1:128  data  rate  expansion  is  first  done  and  is  followed  by 
a 128:1  data  rate  reduction.  The  net  transfer  function  result  should  be  constant  within  the 
passband  constraints  of  the  NT  = 5 design.  A unity-value  sample  followed  by  a string  of 
zeros  was  entered  into  LINT  to  produce  the  impulse  response  whose  frequency  function  is 
that  of  Figure  4-7.  This  impulse  response  data  was  then  entered  into  LDIV.  Note  that  one 
data  point  into  LINT  produces  one  data  point  out  of  LDIV  with  128  points  being  involved  each 
time  at  the  C432/C434  interface. 

This  operation  produced  87  nonzero  final  output  values  which  are  plotted  in  Figure 
5-4.  The  Fourier  transform  of  this  overall  impulse  response  (8192-point  FFT)  is  shown 
in  Figure  5-5.  The  response  is  flat  out  to  f = 0. 9 where  gain  begins  to  drop  to  a -11  dB 
value  (?)  at  f = 1.  Sampling  rate  division  expands  the  frequency  scale  between  input  and 
output.  Note  that  the  input  to  LDIV  had  a band  edge  value  fp  = 0. 00703125  (see  Figure  4-7 
caption);  when  this  f^  value  is  multiplied  by  128  an  fp  = 0. 9 results. 
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Figure  5-5.  Frequency  Response  of  LINT,  LDIV  Chain  for  NT  = J 
(This  is  a Fourier  transform  of  the  Figure  5-4  time  data. 


An  expanded  scale  plot  of  the  Figure  5-5  data  is  shown  in  Figure  5-6 . When  properly 
adjusted  for  normalization  the  "round-trip"  passband  your  variation  limits  are  +0. 0125  dB 
at  f = 0 and  -0. 0226  dB  at  f = 0. 9. 

If  one  were  to  obtain  the  frequency  response  curve  of  the  divider  chain  plotted  in 
terms  of  input  frequency,  the  results  would  be  precisely  those  obtained  by  taking  the 
Fourier  transform  of  the  impulse  response  of  the  equivalent  interpolator  chain.  The  true 
output  frequency  would,  of  course,  be  2L  times  the  input  frequency.  Let  fj  be  the  input 
frequency  and  f^  be  the  divider  output  frequency  (both  normalized  to  R/2).  Define 


f,  . [Vxf,] 

MOD  2 

then 

fd  = ftforftsl 

or 

fd  = 2-ftforft>l 


(5-13) 


(5-14) 


(5-14a) 


The  frequency  division  factors  available  directly  from  C434  are  powers  of  two.  This 

restriction  may  be  lifted  by  use  of  the  interpolation  program  as  a preprocessor  prior  to 

use  of  the  rate  division  program.  Let  the  division  ratio  desired  be  represented  by  which 

is  greater  than  one  but  is  not  a power  of  two.  Let  L^,  Jj  be  the  Interpolation  program 

parameters  as  used  in  Equation  (5-6)  and  related  discussion  in  Section  V4.  par.  5.4.  Let 

L . be  defined  as  the  level  specification  for  the  C434  rate  division  program.  The  division 
a 

ratio  Vd  is  first  specified.  The  Ld,  Lj  and  Jj  values  are  sought  which  will  define  the 
minimum-complexity  interpolator/divider  chain  to  achieve  a division  ratio  Vd,  plus  or 
minus  a tolerance  specification. 


First  choose  the  smallest  Ld  such  that 


(5-15) 


5-15 


If  V'j  is  an  integer,  then 


and 


Li  = Ld 


Ji  - Vd 


/ Integer  \ 
\vd  only  ) 


(5-16) 


which  simply  says  to  multiply  the  sampling  rate  first  by  the  factor  2 a/Vd  in  C432,  then 


divide  this  resulting  sampling  rate  by  2 in  C434. 


If  Vd  is  not  an  integer,  determine  Ld  first  as  in  Equation  (5-15)  then  define  a rate 
multiplication  factor  Kj  of 


and  use  the  procedures  of  Section  V par.  5.4  to  solve  for  L^,  J^,  the  Interpolation  operation 
parameters. 

The  frequency-domain  effect  of  interpolation  is  compression.  A spectrum  that  extends 
from  0 to  (normalized  to  R/2)  is  compressed  to  0 to  Y^/V^  after  a l:Vj  sampling  rate 
increase.  Information  is  not  lost  in  this  process.  The  frequency -domain  effect  of  rate 
reduction  la  expansion  after  filtering.  The  0 to  y1  spectrum  is  first  lowpass  filtered  with 
f'  = y ,/V.  where  V.:  1 is  the  rate  reduction  ratio.  After  this  filtering  operation  the  spec- 

p ' 1 d a 

trum  is  expanded  by  the  factor  Vd  so  that  y^/V^  becomes  y j.  Signal  information  originally 
residing  in  r1/Vd  to  is  lost  in  the  process. 

An  examination  of  the  frequency-domain  effects  of  the  example  of  Equations  (5-15) 
through  (5-17)  will  show  that  the  net  effect  of  Interpolation  followed  rate  reduction  is  an 
equivalent  rate  reduction  operation  of  ratio  Vd  : 1.  If  the  rate  reduction  of  ratio  Vd  precedes 
the  interpolator  of  ratio  Vj  a much  different  result  obtains.  The  original  input  signal  is 
lowpass  filtered  with  f • - yj/Vd.  This  cutoff  frequency  is  then  moved  to  y^  by  the  censoring 
operation.  The  resulting  spectrum  is  then  compressed  by  the  interpolation  ratio  Vj. 
Frequency  components  in  the  original  spectrum  above  y j/Vd  are  suppressed  while  com- 
ponents whose  frequencies  are  below  this  cutoff  value  are  relocated  in  spectral  position  by 
the  factor  V If  Vd  = Vt  the  net  effect  of  this  combined  operation  is  lowpass  filtering 
only.  The  input  and  output  sampling  rates  will  be  equal  and  the  passed  spectral  components 
retain  their  original  frequency  positions.  As  discussed  in  reference  9,  the  above  arrange- 
ment is  an  efficient  approach  to  the  problem  of  lowpass  filtering  when  sampling  rate  changes 
are  not  permitted. 
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SECTION  VI 


COMPUTER  PROGRAMS 

6.1  MULTILEVEL  INTERPOLATION  PROGRAM  LINT 

In  the  Appendix  a computer  program  with  deck  name  C432  is  listed.  Concurrent 
loading  of  the  Block  Data  subprogram  C435  is  required.  The  program  is  written  in 
Honeywell  YFORTRAN  time-share  format.  In  this  section  the  concern  is  with  entry  names 
LINT  and  the  initializing  entry  XLINT,  which  must  be  called  prior  to  the  first  LINT  call. 

CALL  XLINT  (NT,  LEV,  LAG) 

NT  - This  is  an  input  quantity  and  selects  the  Interpolator  type.  Tables  4-1  and  4-2  give 
details.  Types  1-8  assume  significant  spectral  energy  out  to  half  the  sampling  rate, 
and  use  an  extended  y at  the  second  level  in  order  to  maintain  low  stopband  responses. 

NT  values  of  9-12  should  be  used  in  the  oversampled  cases  as  the  y progression  here  is 
a straight  one  to  one-half.  Low  signal  spectrum  levels  are  assumed  in  the  y (R/2)  to  (R/2) 
frequency  range  for  NT  = 9 through  12. 

LEV  - This  input  quantity  is  the  number  of  levels  of  Interpolation.  The  sampling  rate  will 
LEV 

be  multiplied  by  2 . There  is  no  uppler  limit  to  LEV.  However  the  user  must  provide 

LEV 

sufficient  space  in  the  output  array,  F0UT  to  contain  2 values. 

LAG  - This  output  quantity  represents  the  total  delay  through  the  interpolation  processor 
expressed  in  terms  of  samples  at  the  output  sampling  rate. 

After  XLINT  is  used,  a call  to  the  main  entry  LINT  is  made  for  each  input  data  value 

CALL  LINT  (VIN,  F0UT) 

VIN  - This  is  a single  input  sample  value. 

F0UT  - This  is  a vector  in  which  the  output  sample  values  are  returned  by  the  program. 
LEV 

2 values  are  returned  for  each  input  value. 

6.2  MODE  CONVERSION  PROGRAM  C0RE,  REC0.  AND  C0NV 

CALL  C0NV  (NT,  LAGP) 

This  is  an  initiallziiv  entry,  and  must  be  called  before  first  use  of  either  C0RE  or 
REC0.  C0NV  services  both  mode  conversion  routines. 
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NT  - This  is  an  input  quantity  and  selects  the  interpolator  type.  Only  the  first  three  col- 
umns of  Table  4-2  have  any  significance  in  this  usage. 

LAGP  - This  output  quantity  represents  the  mode-conversion  processor  delay  in  terms  of 
output  sample  pairs. 

The  complex-to-real  transformation  is  done  by 
CALL  C0RE  (C0MP,  RE) 

C0MP  - A two-location  input  vector  containing  one  complex  input  data  pair. 

RE  - A two-location  output  vector  containing  a sequential  pair  of  real  samples. 

The  real-to-complex  transformation  is  done  by 
CALL  REC0  (RE,  C0MP) 

RE  - A two-location  input  vector  containing  a sequential  pair  of  real  samples. 

C0MP  - A two-location  output  vector  containing  one  complex  output  data  pair. 

6.3  HP-67  PROGRAM  1-5-79  FOR  SECTION  V.  PAR.  5. 3 EQUATIONS 

This  program  evaluates  Equation  (5-5)  but  in  the  process  it  yields  the  delay  and  N^ 

value  of  Equation  (4-1),  the  Oj  value  of  Equation  (4-2),  the  N^/Oj  ratio  of  Equation  (5-5), 

and  finally  the  N value  of  Equation  (4-3). 
c 

1.  Depress  "D”  to  initialize.  Program  can  handle  level  values  (L)  up  to  15.  All 
{n^}  values  are  set  to  unity  during  this  initialization.  Wait  for  program  halt. 

2.  Enter  number  of  coefficients  at  each  level  starting  with  level  one.  Enter  one 
value  and  depress  "R/S". 

3.  Repeat  step  2 till  {n^}  set  is  stored.  Values  beyond  level  entered  are  one. 

4.  Enter  L value,  depress  "C". 

5.  Program  halts  with  N^  value  of  Equation  (4-1)  in  display.  To  continue 
depress  "R/S". 

6.  Program  halts  with  Oj  of  Equation  (4-2)  in  display.  When  ready  to  continue, 
depress  "R/S". 

7.  Program  halts  with  N^/O^  of  Equation  (5-5)  in  display.  Use  "R/S"  to  continue. 

8.  Program  halts  with  Nc  of  Equation  (4-3)  in  display.  This  completes  the  sequence. 

9.  At  this  point  the  user  can  return  to  step  4 to  repeat  with  a new  L value  and  the 
same  {n^}  set,  or  the  user  may  return  to  step  1 to  enter  new  {n^}  set. 
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6.4  HP-67  PROGRAM  1-11-79  FOR  SECTION  V.  PAR.  5.4  EQUATIQ] 

This  program  may  be  used  to  execute  the  procedure  discussed  in  Section  V,  par.  5.4, 
specifically  with  reference  to  Equations  (5-6),  (5-7),  (5-8),  (5-9). 

1.  Enter  desired  sampling  rate  multiplication  factor  K.  Depress  "A". 

2.  Program  halts  showing  next  L value  to  be  tried.  Press  "R/S"  when  ready 
for  trial  cycle. 

a)  Program  pauses  (flashing  decimal)  and  displays  | A„  | % of  Equation  (5-9). 

b)  Program  pauses  to  display  J value  of  Equation  (5-7). 

A 

c)  Program  pauses  to  display  K value  of  Equation  (5-8). 

3.  Program  increments  L and  returns  to  step  2. 

Note:  During  step  2,  halt,  user  may  enter  any  L value  for  trial.  If  user 

enters  L value  and  "STC^'  L sequence  is  modified.  For  data  recovery, 
incremented  L is  in  register  2,  last  J value  is  in  register  1. 

6.5  MULTILEVEL  SAMPLING  RATE  DIVISION  PROGRAM  LDIV 

This  program  is  listed  under  deck  name  C434  and  requires  concurrent  loading  of  the 
Block  Data  subprogram  C435.  The  program  is  written  in  Honeywell  YFORTRAN  time-share 
format.  There  are  two  entry  points  to  this  program.  An  initializing  entry  XLDTV  must  be 
called  prior  to  the  first  usage  of  the  working  entry  LDIV. 

CALL  XUMV  (NT,  LEV,  LAG) 

NT  - An  input  variable  which  selects  the  processor  type.  Type  specifications  1 through  12 
are  currently  valid. 

LEV  - An  input  variable  which  specifies  the  number  of  levels  of  rate  division.  The 

I FV 

sampling  rate  will  be  divided  by  2 . There  is  no  upper  limit  to  LEV,  however  the 

LEV 

user  must  dimension  an  input  array  FIN  of  size  2 or  greater. 

LAG  - This  output  variable  represents  the  total  delay  through  the  rate  division  processor 
expressed  in  terms  of  input  samples. 

After  XLDIV  is  called,  a call  to  the  main  entry  LDIV  is  made  for  each  group  of 

LEV 

2 input  values;  one  output  value  results  for  each  LDIV  call. 

CALL  LDIV  (FIN,  VOUT) 

LEV 

FIN  - An  input  vector  containing  2 data  values. 

VOUT  - A single  output  variable. 
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CONCLUSIONS 

The  ease  of  design,  computational  economy,  and  high  data-rate  increase  factors 
available  from  the  iterative-interpolation  technique  make  this  approach  attractive  relative 
to  the  conventional  filter  method  of  sample-rate  multiplication.  The  filter  synthesis 
procedure  which  was  derived  from  the  interpolation  algorithm  also  appears  to  have 
practical  utility. 

The  computer  program  C432  of  the  Appendix  should  be  useful  in  its  present  form  for  a 
variety  of  interpolation,  filter  design,  and  mode-conversion  function.  Some  improvements 
in  the  program  are  planned,  however.  A y = 95%  cascade  should  be  added.  This  gives  a 
transition  region  which  is  11%  of  the  passband.  Such  designs  are  unavoidably  expensive 
in  processor  length.  If  digital  filters  are  ever  to  rival  analog  filters  in  performance,  hard- 
ware technology  must  advance  to  the  point  where  such  designs  can  be  accepted.  The  rate- 
doubling sections  of  the  cascades  should  be  increased  in  number  to  accommodate  some  lower 
stopband  level  specifications.  As  previously  discussed,  the  linear  interpolation  stage  should 
be  recoded  for  an  arbitrary  1:J  sampling  rate  increase.  The  algorithms  in  Section  V 
par.  5.4  for  arbitrary  sampling  rate  multiplication  factors  could  be  included  in  a future 
version  of  the  program.  The  user  would  then  be  required  only  to  specify  the  sampling 
rate  multiplication  factor  desired  and  the  program  would  automatically  provide  the  optimum 
solution. 

Notice  should  be  taken  of  the  real  advantages  of  applying  the  iterative  interpolation 
technique  as  opposed  to  the  more  conventional  use  of  low-pass  or  band  stop  filters  for  this 
service.  The  ease  with  which  good  Interpolators  can  be  designed  and  their  efficiency  in  ser- 
vice, especially  at  high  multiplication  ratios,  should  suggest  many  new  areas  of  application 
for  the  interpolation  process. 

One  such  new  area  of  promise  concerns  the  use  we  have  made  here  of  the  Interpolator 
in  filter  design.  Basically,  a generic  impulse  response  function,  believed  to  be  optimum 
in  some  sense,  primes  the  process.  Whether  this  generic  function  is  obtained  out  of  first 
level  as  in  our  example,  is  or  supplied  externally  is  incidental.  The  proper  interpolation 
process  is  then  applied  to  derive  filters  having  passband  frequency  functions  which  are 
scaled  replicas  of  the  generic  filter  equivalent  functions.  It  would  be  useful  to  know  to  what 
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degree  the  original  optimality  property  is  retained  in  the  derived  filters.  If  the  degree  of 
retention  is  high,  these  results  could  have  significance  beyond  the  digital  filter  design 
goals  of  Section  IV. 


The  complementary  use  of  the  cascade  technique  for  sampling  rate  division  In  pro- 
gram C434  makes  available  many  of  the  advantages  of  the  companion  C432  Interpolation 
program.  A combination  of  these  two  programs  offers  the  user  a convenient,  efficient  set 
of  procedures  for  a variety  of  signal  processing  tasks. 
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APPENDIX  A 


PROGRAMS 

In  this  Appendix  four  computer  program  listings  are  given  as  mentioned  in  the 
main  text.  Block  data  subprogram  C435  must  be  loaded  along  with  C432  or  C434.  The 
two  HP-67  programs  have  been  coded  so  that  they  may  be  merged  onto  one  magnetic  card, 
if  desired.  The  FORTRAN  programs  are  also  available  in  standard  (FORTY)  source-deck 
format. 


A-l 


10*C432  C432  ITERATIVE  I NTERPOL ATI UN  SUBROUTINE 

20*  JP  COSTAS  12/ 13/ 18 

30*USE  WITH  BLOCK  DATA  SUBPROGRAM  C435 
40* 

50  SUBROUTINE  LINT( VI N# POUT) 

40* 

70*VIN-SINGLE  INPUT  DATA  ITEM 
80*F0UT-ARRAY  FOR  OUTPUT  DATA  STORAGE 
90* 

100* DIMENSIONED  FOR  MAX  N OF  44 
1 10* 

120  COMMON  /GP432/  SCOEFC  512)  #NGC  6#20) 

130* 

140  DIMENSION  F0UT<  t)#H0LDC  128#5)#0UT<  32#  6>#NCF<  S)#MP<  5)# 
I 504  C0EFC  44#  5)  # DELAY ( 1 28  ) # USE*  128)#  C0EF2C  64)  # 

1604  COMPC  2) # RE<  2) 

1 70* 

180  E0U1  VALENCE  ( DELAY < 1 > #HOLD(  1 • 1 ) >#  ( USE*  1 )#HOLD(  l#2)  ) • 
1904  ( C0EF2C  1 ) # C0EF<  1 # 1 ) ) 

200* 

210* 

220*NT- INTERPOLATION  TYPE  1-12  IN  THIS  VERSION 
230*NGC6#NT) -ADDRESS- 1 IN  SCOEF  OF  COEFFICIENT  DATA 
240*NG(L#NT) -NUMBER  OF  COEFS  AT  LEVEL  L 
250* 

260  0UT<1#1)»V1N 

270  DO  1 L«1#NLEV 
280  NCFL-NCFCL) 

290  NCFM«NCFL-1 
300  LENL«2*NCFCL) 

310  LML«LENL-1 
320  KMAX*2**<L-1 ) 

330*KMAX  IS  NUM  INPUT  DATA  VALUES  AT  THIS  LEVEL 

340  KK»0 

350  DO  1 K* 1 # KMAX 

360  MP(L) *M0D< MP<  L) ♦ 1 #LENL) 

370  MPL-MPCL) 

380  KK«KK*1 

390  NEW«MOD<MPL*NCFM#LENL) 

400  HOLDC  NEW*  1 #L>  *0UT<  K#L) 

410  B0UT*O* 

420  M1»MPL 

430  M2«M0D<  MPL ♦LML • L ENL ) 

440  DO  2 J-WNCFL 
450  M1*M0D<M1*LML#LENL) 

460  M2>M0D<N2*1#LENL> 

470  2 B0UT*B0UT*COEF<  J# L> *<H0LD< Ml  ♦ 1 #L)  *H0LD< M2*  1 #L)  ) 


I 


I 


480  0UTC  KK*L+1 )*B0UT 
490  KK«KK*1 

500  1 0UT<  KK#L* 1 ) *H0LD< MPL* 1 »L) 

510  1 F<LEV* 6T*  5)  GO  TO  3 

520  DO  4 K»  1 * KA 

530  4 FOUT<K)>0UT(K*NA) 

540  RETURN 

550  3 KK-0 

540  DO  4 K»1 * 32 

570  B0UT«<0UTCK*6)-SAVE)/DEL 

580  DO  5 J*1«NDEL 

590  KK*KK*1 

400  5 FOUTt  KK» *SAVE+FL0AT( J~1 ) 6B0UT  . 

610  6 SAVE*OUT<  K*  6) 

620  RETURN 

6306 

6404 

650  ENTRY  XLINTCNT*LEV*LAG) 

6604 

6704 

6804IN1 T1 AL1 ZAT10N  ENTRY  FOR  LINT 

4906NT-INTERP0LATI0N  TYPE  (1-12  IN  THIS  VERSION) 

7004LEV*NUM  LEVELS  OF  INTERPOLATION 

7 10* EACH  INPUT  VALUE  GENERATES  244LEV  OUTPUT  VALUES 

7204L A G~ DELAY  THRU  INTERPOLATOR  IN  TERMS  OF  OUTPUT  SAMPLES 

7306 

7406 

750  SAVE«0. 

740  KK«NGC6*NT> 

770  DO  7 L*l*  5 
780  MPCD--1 
790  M1«NG<L#NT) 

800  NCF(L)«M1 
810  DO  7 J* 1 * Ml 
820  KK*KK*1 

830  7 C0EFCJ*L)«SC0EF(KK)/2. 

8 40  NLEV*M1N0<LEV*5) 

8 50  DO  8 J«  1*440 
840  8 HOLDC  J*  1>*0* 

8 70  DO  9 J*l*  192 
880  9 0UT< J* 1 )*0» 

890  I FCLEV*  6T»5)  GO  TO  10 

900  KA*2*4NLEV 

910  NA*NLEV*1 

920  GO  T01 4 

930  10  NDEL*266CLEV*5) 

9 40  DEL-FLOATiNDEL) 


A-3 


950  1 * Ml  ■LEV*! 

9 40  LAG>1-8**LEV 
9 70  DO  15  J-WLEV 
980  KK*l 

990  IF<  J«LE*5)  KK*NG< J*NT» 

1000  15  LA6-LAG*KK*8**<M1-J> 

1010  RETURN 
1080* 

1030* 

1040  ENTRY  CONVCN  T»LAGP> 

1 050*1 NI TI AL1ZIN6  ROUTINE  FOR  BOTH  CORE  AND  RECO 
1060*NT*INTERP8LAT10N  TYRE  1*18  IN  THIS  VERSION 
1 070*LAGP 'PROCESSOR  DELAY  IN  TERMS  OF  OUTPUT  SAMPLE  PAIRS 


108b  NCPL-NGC  l«NT> 

1090  LA6P-NCFL 
1100  LENL>8*NCFL 
1110  LML-LENL-1 
1180  MPL— 1 
1130  KK»N6C  &»NT> 

1140  DO  11  J»|«NCFL 
1 • 50  KK«KK*I 

11*0  11  C0EF8( J) ■SCOEFC  KK) /8» 

1170  DO  18  J«W85* 

1180  18  HOLDC J#  1>*0» 

1190  FLIP-1. 

1 800  RETURN 
1810* 

1880* 

1830  ENTRY  COREC  COMP#  RE> 

1 840*C0MPLEX -TO'REAL  TRANSFORMATION  ENTRY 
1 S50*C0MP-L8CATI ON  PAIR  FOR  COMPLEX  INPUT 
1 8*O*RE'L0CATI0N  PAIR  F0R  REAL  OUTPUT 
1870  MPL *M0 D(  NPL ♦ I * LENL > 

1880  RE< I >«DELAY<MPL*1 >*FLIP 

1890  FLIP-FLIP 

1300  NEW>M0D(MPL*NCFL*LENL> 

1310  DELAYS  NEW* 1 > *COMP(  1 ) 

1380  USECNEW+l )*C0MPC  8> 

1330  B0UT-O* 

1340  M1«M0D(MPL*I»LENL> 

1 350  M8-MPL 

1360  D0  13  J«1*NCFL 

1370  Ml aM0D<  Ml *LML»LENL> 

1380  M2*M0D<  M2*  1 # LENL) 

1390  1 3 B0UTaB0UT*C0EF2< J)*<  USE(M| *1 ) *USE(M8*I > > 
1400  RE<2>  "BOUT*  FLIP 

1410  RETURN 
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UNCLASSIFIED 


GENERAL  ELECTRIC  CO  SYRACUSE  N Y HEAVY  MILITARY  EGUI—ETC  F/6  9/5 
AN  ITERATIVE  INTERPOLATION  TECHNIQUE  USING  FINITE  IMPULSE  RESPO— ETC(U) 
FEB  79  J P COSTAS 

R79EMH1  NL 


1 


1480* 

1430* 

1440  ENTRY  REC0<RE»C0MP) 

I 4S0*RE -LOCATION  PAIR  FOR  REAL  INPUT 
1 440*C0NP -LOCATION  PAIR  FOR  COMPLEX  OUTPUT 
1470  HPL»MO!XMPL*t»LENL> 

1400  C0NP<8>«  DELAY  <MPL*I> 

1490  NE  W*MOCK  MPL+NCFL#LENL) 

1500  USE<NEW*1>«RE<1>*FLIP 
1510  FLIP— FLIP 
1580  DELAY C NEW* 1>«RE<8>* FLIP 
1530  B0UT>0» 

1540  H1«M0D<MPL*ULENL> 

1 550  N8-MPL 
1 560  DO  14  J«1#NCFL 
.1570  M1»M0D<M1*LML»LENL) 

I 500  M8«N0D<  M8*  1 * LEND 

1590  14  B0UT*B0UT*C0EF8<  J> •<  USEC N 1 * I > *USE<  N8* 1 > > 
1400  COMPC 1 > »BOUT 
1410  RETURN 
I 480* 

1430* 

1440  END 
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104C434  C434  ITERATIVE  SAMPLING  RATE  DIVISION  SUBROUTINE 

064JP  COSTAS  2/1/19 

304 USE  WITH  BLOCK  DATA  SUBPROGRAM  C435 
404 

SO  SUBROUTINE  LDI  VC  FIN#  VOUT) 

004 

70*FIN- INPUT  DATA  ARRAY 
804V0UT-S1NGLE  OUTPUT  DATA  VALUE 
904 

10O4DIMENS1ONED  FOR  MAX  N OP  04 
1104 

100  COMMON  /GP430/  SCOEPC  510>«NG<  4»00> 

1304 

140  DIMENSION  FIN<  1>#  DELAY ( 6S«S>»H0LD(  1O8»S)«0UT< 30#  4)#NCP(5>* 
1500  MPC  5)*MPD<  5)»C0EF<  04*5) 

1004 

I 704LAST  EFN  IS  17 

1004 

1904 

OOO4NT-1NTERPOLAT10N  TYPE  1-10  IN  THIS  VERSION 
0104NGC  0# NT) "ADDRESS- 1 IN  SCOEF  OF  COEFFICIENT  DATA 
0004NGC  L# NT) "NUMBER  OF  COEPS  AT  LEVEL  L 
0304 

040  IFfLEV.LC*  5)  60  TO  3 

0504 ENTER  LINEAR  INTERPOLATION  PROCESSOR 

000  SUN-FINC I ) 

070  DO  0 J«1#NDELM 

000  0 SUM«SUM*FLOATC NOEL -4) /DEL4F1NC J*  I > 

090  0UT< 1# 1)>SUM*SAVE 
300  MPL-1. 

310  DO  II  K«0»30 
300  MPL-MPL+NDEL 
330  SUM-FINC  MPL) 

340  DO  10  <J*!«NDELM 

•350  10  SUM«SUM^FL0AT<NDEL-J>/DEL4CFINCMPL«J>^FIN(MPL«J>> 

300  11  OUTCK* 1)*SUM 

3704N0H  SET  UP  S*VE  FOR  NEXT  CYCLE 

380  SAVE-O. 

390  DO  13  J»1#NDELM 

400  13  SAVE«SAVE*FL0ATCJ>/DEL4F1N<KB*J> 

410  GO  TO  14 

400  3 DO  15  K-I«KA 

430  1 5 OUTCK*  1)-FIN<K> 

440  14  CONTINUE 

4504 

4404 

4 704 ENTER  THE  NORMAL  CASCADE  PROCESSOR 
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480 

DO  1 L-1»NLEV 

490 

LP-NLEV*1-L 

500 

KMAX«8**LP 

510*KHAX  IS  NUM  INPUT  DATA  VALUES  AT  THIS  LEVEL 

seo 

NCPL>NCP(LP) 

530 

LENL>8*NCPL 

540 

LEND-NCFL*! 

550 

LHL-LENL-I 

5*0 

KK-0 

5704KK  IS  POINTER  IN  DATA  OUTPUT  ARRAY 

580* 

590 

DO  1 K-1»KNAX»8 

400 

HP<L)-H0D<HP<L)+1»LENL) 

410 

NPL-HPCL) 

480 

MPD(L)-M0D(MPD(L)«1»LEND) 

430 

MPDL-NPD(L) 

440 

KK-KK*I 

450 

NEH"NOD( NPL*NCPL*LENL) 

440 

NEVD*NOD<  MPDL^NCFL#  LEND) 

470 

HOLD(NE*H#L)-0UT<K«L) 

480 

DELAY(NEWD* 1 #L> «OUT< K* 1 #L) 

490 

B0UT»DELAY( HPCL* 1 • L ) 

700 

H1-H0D(HPL+I«LENL) 

710 

M8»MPL 

780 

DO  8 J- 1 »NCFL 

730 

M 1 »M0D< N 1 ♦LNL* LENL ) 

740 

M2-M0DC  M2*  1 » LEND 

750  i 

2 B0UT-B0UT*C0EP<J*LP)*<H0LD<M1*|»L)*H0LD<M8*1»L)) 

740 

1 OUT(KK»L*1)«B0UT 

770 

VOUT-BOUT, 

780 

RETURN 

790* 

800* 

810 

ENTRY  XLDI V(NT#LEV#LA6) 

880* 

830* 

8 40*1  NITIALI ZING  ENTRY  FOR  LDIV 

8 S04NT-PR0CESS8R  TYPE  <1-18  IN  THIS  VERSI 8N>  (INPUT  DATA! 

840*LEV-NUM  LEVELS  8P  DIVISION  (INPUT  DATA) 

8 70*LA6-PR0CESS0R  DELAY  IN  TERNS  OP  INPUT  PULSES  (OUTPUT  DATA) 

880*0NE  OUTPUT  VALUE  RE0U1RES  8**LEV  INPUT  VALUES 

890* 

900 

SAVE-O* 

910 

KK-NG(  4#NT) 

980 

DO  7 L-l#5 

930 

HPCL)  — 1 

940 

NPD(L)  — 1 

990  Mt»NGCL»NT> 

940  NCFCD-Ml 
970  DO  7 
980  KK-KKM 

990  7 COCP(J»L>«SCOtr(KK)/e. 
1000  NLEV*MIN0!LEV#9> 

1010  DO  t J*l# 440 
1020  8 HOLD!  J#  1>*0» 

1030  DO  9 J-1.192 
1040  9 OUTC J# 1)*0« 

1090  DO  14  J-W329 
1040  14  DELAY! 

1070  NDEL*244!  LEV-5) 

1080  DEL* FLOAT! N DEL > 

1090  KA«£**LEV 
1100  KB*KA-NDCL*1 
1110  NDELH-NDFL-1 
1120  M1*LEV«1 
1130  LA 6* 1 “244LEV 
1140  DO  17  J*1»LEV 
1190  KK*| 

1140  IP!J*LE*9)  KK*H6! J#NT> 
1170  17  LA6*LAG4KK«2*«!N1-J) 
1 180  RETURN 
11904 
1200  END 


480 8.  I • 820 58  . - • 88  5339  . . 0 75*  446. 

49081  • 1 4551  *-•  1 48007# 

50081  *01  401 . 

51081  .00348/ 

580* 

530*NT*4.  80*.  -85DB 

540  DATAC SCOCPC  J>.«JaS1.58>/l  *2588  6. -*366063.  *814951# 

550* 

5 8081  • 8085  7. -•  854078. 

57081*05784. 

58081*01401. 

59081*00348/ 

800* 

8I0*NT*5.  90*.  -65DB 
620 DATAC  SC0EFCJ3.  J-59.89) /I  *87098. -• 41  7871.  *843548. -•  1 86887.  . 188307. 
630  4-9 • 300 66E-08. 7*199  47E-08. -5 • 61 48 1 E-08. 4*  38098E-02. -3 • 408 1 5E-08. 
64088*  818  41 E-08. - 1 *98938E-08. I • 48  784E-08. - 1 *088  60E-08. 7*  7488  4E-03. 
6504-S*  3 449E-03.  3*  S1893E-03.  -2*  I8409E-03.  1 • 70604E-03. 

680* 

67081 *83918. -*331 141. *184558. -*040 7535.8* 750 48E-03. 

680  81  • 18469. -*815619.  3*  1049  6E-08. 

69041 *13018. -*130389. 

70081  *18688. -*18689/ 

710* 

720*NT»6.  90S.  -50DB 

7 30 DAT A<  SCOEFC  J)  # J»90.  1 18J/1  *27034. -.41577*  *840448.-*!  68436.  *1 1 71  66. 
7408-8  • 70883E-02.  6* 54880E-02.  -4.92510E-02.  3*  67355C-08.  -2.69847E-02. 
75081 *93  773E-08. - 1 *3457  6E-02.8  *9  4S8E-03.  -8*51 445E-03. 

760* 

77081 *83183. -*31873. *103554. -*0851 79. 

78081  • 1 4551.  -•  1 48007. 

79081 *13018. -.130 389. 

80081*00348/ 

810* 

880*NT"7.  90S.  -40DB 

830DATA<SC0EP<J>.J*I  13.189) /l  *86955. -.41 3388.  *83656. -.1571 44.  . 

8 404* 1 10789.  -.07976!  5.  *0576061  • -*04101 65.  *088  7088. -*031883. 

850* 

8 6041  *88058. -*88  5339  • .0759  446. 

8 7041  * 1 4551. -*  1 48007. 

88081*01401. 

89081*00348/ 

900*NT-8.  90S.  -8508 

9 1 OOATAC  SCOEFC  J)  . 1 30. 1 40 > /!  * 268  44.  - * 409  63  7.  . 8309  61  • - • 1 49  668. 

9804*  1031 41. -.1 31 136. 

930* 

94061  *80857.  -*854076. 
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9 50*1*05764# 

960*1*01401# 

970*1  *00348/ 

900* 

990*NT*9#  (OX*  -6508 

1 QOODATAC  SCOCrC  J)  » J«  1 41  # 1 53>  / 1 • 239 18  • - . 33 1 1 41  # * 1 84558  • - • 040  7535# 
1 010*. 008  75048# 

1020* 

1030*1  •184*9#—  815*19#  .031049*# 

1040*1 *13018* ••130389# 

1 050*1  • 1 8*28# -•  1 2609# 

1060*1  .00087/ 

1070* 

1080*NT«10#  60S*  *5008 

1090DATA< 3C0tr<J>#J«l 54# 163>/l *83183# -*31873# *103554# -*0851 79# 

1 100*1 .14551#-#  1 48007# 

11 10*1 *13018# -#130389# 

1180*1*00348# 

1130*1*0008 7/ 

1148* 

I !50*NT-1 I#  608#  -4008 

1 1 60DATAC  SCiCFC  J)  # J«1 64#  1 71 ) /l  • 88058 • - *88 5339#  *0  759  446# 

1170* 

1 180*1  *14551  #-*148007# 

1190*1*01401# 

1800*1*00348# 

1810*1*00087/ 

1880* 

1 830*NT* 1 8#  60S*  -8508 

1 2 40 OAT* < SCOCPC  J 1 # J«  1 78#  1 7 7>  / 1 . 8085 7#  - * 8540  76# 

1850* 

1860*1  *05764# 

1870*1*01401# 

1880*1*00345# 

1890*1*00067/ 


PROGRAM  1-5-79  (HP-67),  FOR  USE  WITH  SECTION  V,  PAR.  5.3  EQUATIONS 


1 

<D) 

FLBLD 

2 

2 

3 

5 

4 

hSTI 

5 

(5) 

fLBL5 

6 

£DSZ 

7 

1 

8 

STO(i) 

9 

hRCI 

10 

1 

11 

0 

12 

gx/y 

13 

GT05 

14 

(3) 

fLBL3 

15 

R/S 

Enter  (N.) 

K 

16 

STO(l) 

17 

asz 

18 

GTOS 

19 

(C) 

fLBLC 

Enter  L 

20 

STOO 

21 

2 

22 

hx  ~y 

23 

by*  t 

24 

ENTt 

25 

ENTt 

26 

1 

27 

- 

28 

ST02 

29 

hRi 

30 

• 

2 

31 

X 

32 

ST03 

33 

CLX 

34 

ST05 

35 

ST  06 

36 

• 

37 

5 

38 

ST  04 

39 

9 

40 

hSTI 

41 

(2) 

ILBL2 

42 

flSZ 

43 

RCL4 

44 

2 

45 

X 

46 

ST  04 

iV 

47 

48 

RCL(i) 

X 

49 

ST  0+6 

0.  sum 

SO 

RCL3 

1 

51 

2 

52 

-r 

53 

STOS 

54 

RCL(i) 

55 

X 

56 

STO+5 

N-  sum 

57 

hRCI 

I 

58 

9 

59 

- 

1-9 

60 

RCLO 

61 

gx/ty 

62 

GTO2 

63 

RCL5 

64 

RCL2 

65 

- 

66 

R/S 

Show  N«  D. 

67 

RCL6 

I 1 

68 

R/S 

Show  0. 

69 

• 

1 

70 

R/S 

Show  N^/Oj 

71 

RCL5 

72 

2 

73 

RCLO 

74 

1 

75 

_ 

76 

h/ 

77 

■p 

78 

ENTt 

79 

QNT 

80 

ST07 

81 

gx=y? 

82 

GT04 

83 

1 

84 

8TO+7 

85 

(4)  tLBLA 

86 

RCL7 

87 

R/S 

Show  N 

0 
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